Changeset 5772


Ignore:
Timestamp:
09/13/10 13:18:07 (15 years ago)
Author:
Eric.Larour
Message:

New solution strategy, based on creation of Kff, Kfs directly,
instead of Kgg. Idem for pf instead of pg.
Kept backwards compatibility for now, and only applied to macayeal diagnostic
in 2d. Will switch to this new strategy in the near future.

Location:
issm/trunk/src
Files:
13 added
48 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/Constraints.cpp

    r5057 r5772  
    6363}
    6464/*}}}*/
    65 /*FUNCTION Constraints::SetupSpcs{{{1*/
    66 void   Constraints::SetupSpcs(Nodes* nodes,Vec yg,int analysis_type){
    67 
    68         int i;
    69 
    70         Node* node=NULL;
    71         int nodeid;
    72         int dof;
    73         double value;
    74 
    75         for(i=0;i<this->Size();i++){
    76 
    77                 Object* object=(Object*)this->GetObjectByOffset(i);
    78 
    79                 /*Check this is a single point constraint (spc): */
    80                 if(object->Enum()==SpcEnum){
    81 
    82                         Spc* spc=(Spc*)object;
    83 
    84                         if(spc->InAnalysis(analysis_type)){
    85 
    86                                 /*Ok, this object is a constraint. Get the nodeid from the node it applies to: */
    87                                 nodeid=spc->GetNodeId();
    88                                 dof=spc->GetDof();
    89                                 value=spc->GetValue();
    90 
    91                                 /*Now, chase through nodes and find the corect node: */
    92                                 node=(Node*)nodes->GetObjectById(NULL,nodeid);
    93 
    94                                 /*Apply constraint: */
    95                                 if(node){ //in case the spc is dealing with a node on another cpu
    96                                         node->ApplyConstraint(yg,dof,value);
    97                                 }
    98                         }
    99 
    100                 }
    101         }
    102 
    103         /*Assemble yg: */
    104         VecAssemblyBegin(yg);
    105         VecAssemblyEnd(yg);
    106 }
    107 /*}}}*/

  • issm/trunk/src/c/Container/Constraints.h

    r5057 r5772  
    2828                /*numerics: {{{1*/
    2929                int   NumberOfConstraints(void);
    30                 void  SetupSpcs(Nodes* nodes,Vec yg,int analysis_type);
    3130                /*}}}*/
    3231

  • issm/trunk/src/c/Container/Nodes.cpp

    r5296 r5772  
    4444/*Numerics*/
    4545/*FUNCTION Nodes::DistributeDofs{{{1*/
    46 void  Nodes::DistributeDofs(int analysis_type){
     46void  Nodes::DistributeDofs(int analysis_type,int setenum){
    4747
    4848        extern int num_procs;
    … …  
    5858        int  numnodes=0;
    5959
     60        /*some check: */
     61        if ((setenum!=GsetEnum) && (setenum!=FsetEnum) && (setenum!=SsetEnum))ISSMERROR("%s%s%s"," dof distribution for set of enum type ",EnumToString(setenum)," not supported yet!");
     62
    6063        /*Go through objects, and distribute dofs locally, from 0 to numberofdofs: */
    6164        for (i=0;i<this->Size();i++){
    … …  
    6467                /*Check that this node corresponds to our analysis currently being carried out: */
    6568                if (node->InAnalysis(analysis_type)){
    66                         node->DistributeDofs(&dofcount);
    67                 }
    68         }
    69 
    70 
     69                        node->DistributeDofs(&dofcount,setenum);
     70                }
     71        }
    7172
    7273        /*Ok, now every object has distributed dofs, but locally, and with a dof count starting from
    … …  
    9697                Node* node=(Node*)this->GetObjectByOffset(i);
    9798                if (node->InAnalysis(analysis_type)){
    98                         node->OffsetDofs(dofcount);
     99                        node->OffsetDofs(dofcount,setenum);
    99100                }
    100101
    … …  
    104105         * object that is not a clone, tell them to show their dofs, so that later on, they can get picked
    105106         * up by their clones: */
    106         maxdofspernode=this->MaxNumDofs(analysis_type);
     107        maxdofspernode=this->MaxNumDofs(analysis_type,setenum);
    107108        numnodes=this->NumberOfNodes(analysis_type);
    108109
    … …  
    114115                Node* node=(Node*)this->GetObjectByOffset(i);
    115116                if (node->InAnalysis(analysis_type)){
    116                         node->ShowTrueDofs(truedofs,maxdofspernode);//give maxdofspernode, column size, so that nodes can index into truedofs
     117                        node->ShowTrueDofs(truedofs,maxdofspernode,setenum);//give maxdofspernode, column size, so that nodes can index into truedofs
    117118                }
    118119        }
    … …  
    125126                Node* node=(Node*)this->GetObjectByOffset(i);
    126127                if (node->InAnalysis(analysis_type)){
    127                         node->UpdateCloneDofs(alltruedofs,maxdofspernode); //give maxdofspernode, column size, so that nodes can index into alltruedofs
     128                        node->UpdateCloneDofs(alltruedofs,maxdofspernode,setenum); //give maxdofspernode, column size, so that nodes can index into alltruedofs
    128129                }
    129130        }
    … …  
    215216/*}}}*/
    216217/*FUNCTION Nodes::NumberOfDofs{{{1*/
    217 int   Nodes::NumberOfDofs(int analysis_type){
     218int   Nodes::NumberOfDofs(int analysis_type,int setenum){
    218219
    219220        int i;
    … …  
    233234                        if (!node->IsClone()){
    234235
    235                                 numdofs+=node->GetNumberOfDofs();
     236                                numdofs+=node->GetNumberOfDofs(NoneApproximationEnum,setenum);
    236237
    237238                        }
    … …  
    241242        /*Gather from all cpus: */
    242243        MPI_Allreduce ( (void*)&numdofs,(void*)&allnumdofs,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
    243 
    244244        return allnumdofs;
    245245}
    … …  
    335335/*}}}*/
    336336/*FUNCTION Nodes::MaxNumDofs{{{1*/
    337 int   Nodes::MaxNumDofs(int analysis_type){
     337int   Nodes::MaxNumDofs(int analysis_type,int setenum){
    338338
    339339        int i;
    … …  
    350350                if (node->InAnalysis(analysis_type)){
    351351
    352                         numdofs=node->GetNumberOfDofs();
     352                        numdofs=node->GetNumberOfDofs(NoneApproximationEnum,setenum);
    353353                        if (numdofs>max)max=numdofs;
    354354                }

  • issm/trunk/src/c/Container/Nodes.h

    r5114 r5772  
    1919                /*}}}*/
    2020                /*numerics: {{{1*/
    21                 void  DistributeDofs(int analysis_type);
     21                void  DistributeDofs(int analysis_type,int SETENUM);
    2222                void  FlagClones(int analysis_type);
    2323                void  FlagNodeSets(Vec pv_g, Vec pv_f, Vec pv_s,int analysis_type);
    24                 int   NumberOfDofs(int analysis_type);
     24                int   NumberOfDofs(int analysis_type,int setenum);
    2525                int   NumberOfNodes(int analysis_type);
    2626                int   NumberOfNodes(void);
    2727                void  Ranks(int* ranks,int analysis_type);
    28                 int   MaxNumDofs(int analysis_type);
     28                int   MaxNumDofs(int analysis_type,int setenum);
    2929                /*}}}*/
    3030

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

    r5740 r5772  
    172172                                        ./objects/Loads/Numericalflux.cpp\
    173173                                        ./objects/Loads/Numericalflux.h\
     174                                        ./objects/Numerics/ElementMatrix.h\
     175                                        ./objects/Numerics/ElementMatrix.cpp\
     176                                        ./objects/Numerics/ElementVector.h\
     177                                        ./objects/Numerics/ElementVector.cpp\
    174178                                        ./objects/Params/Param.h\
    175179                                        ./objects/Params/BoolParam.cpp\
    … …  
    448452                                        ./modules/CostFunctionx/CostFunctionx.h\
    449453                                        ./modules/CostFunctionx/CostFunctionx.cpp\
     454                                        ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h\
     455                                        ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp\
    450456                                        ./modules/Orthx/Orthx.h\
    451457                                        ./modules/Orthx/Orthx.cpp\
    … …  
    524530                                        ./modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h\
    525531                                        ./modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\
     532                                        ./modules/Reduceloadx/Reduceloadx.h\
     533                                        ./modules/Reduceloadx/Reduceloadx.cpp\
    526534                                        ./modules/NodeConnectivityx/NodeConnectivityx.cpp\
    527535                                        ./modules/NodeConnectivityx/NodeConnectivityx.h\
    … …  
    723731                                        ./objects/Loads/Numericalflux.cpp\
    724732                                        ./objects/Loads/Numericalflux.h\
     733                                        ./objects/Numerics/ElementMatrix.h\
     734                                        ./objects/Numerics/ElementMatrix.cpp\
     735                                        ./objects/Numerics/ElementVector.h\
     736                                        ./objects/Numerics/ElementVector.cpp\
    725737                                        ./objects/Params/Param.h\
    726738                                        ./objects/Params/BoolParam.cpp\
    … …  
    9961008                                        ./modules/CostFunctionx/CostFunctionx.h\
    9971009                                        ./modules/CostFunctionx/CostFunctionx.cpp\
     1010                                        ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h\
     1011                                        ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp\
    9981012                                        ./modules/Orthx/Orthx.h\
    9991013                                        ./modules/Orthx/Orthx.cpp\
    … …  
    10701084                                        ./modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h\
    10711085                                        ./modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\
     1086                                        ./modules/Reduceloadx/Reduceloadx.h\
     1087                                        ./modules/Reduceloadx/Reduceloadx.cpp\
    10721088                                        ./modules/MassFluxx/MassFluxx.cpp\
    10731089                                        ./modules/MassFluxx/MassFluxx.h\
    … …  
    11281144                                        ./solvers/solver_adjoint_linear.cpp\
    11291145                                        ./solvers/solver_diagnostic_nonlinear.cpp\
     1146                                        ./solvers/solver_diagnostic_nonlinear_kgg.cpp\
     1147                                        ./solvers/solver_diagnostic_nonlinear_kff.cpp\
    11301148                                        ./solvers/solver_thermal_nonlinear.cpp\
    11311149                                        ./modules/Bamgx/Bamgx.cpp\

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

    r5057 r5772  
    2929
    3030        /*Get gsize; */
    31         gsize=nodes->NumberOfDofs(analysis_type);
     31        gsize=nodes->NumberOfDofs(analysis_type,GsetEnum);
    3232
    3333        if(gsize){

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

    r5103 r5772  
    2525
    2626        /*Get size of vector: */
    27         gsize=nodes->NumberOfDofs(configuration_type);
     27        gsize=nodes->NumberOfDofs(configuration_type,GsetEnum);
    2828        if (gsize==0) ISSMERROR("Allocating a Vec of size 0 as gsize=0 for configuration: %s",EnumToString(configuration_type));
    2929       

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

    r5578 r5772  
    6666        parameters->AddObject(new StringParam(SolverStringEnum,iomodel->solverstring));
    6767        parameters->AddObject(new IntParam(NumberOfElementsEnum,iomodel->numberofelements));
     68        parameters->AddObject(new BoolParam(KffEnum,iomodel->kff));
    6869
    6970        /*Deal with more complex parameters*/

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

    r5114 r5772  
    2626                 *a  node has already been distributed dofs on one cpu, all other cpus with the same node cannot distribute it
    2727                 *anymore. Use clone field to be sure of that: */
    28                 nodes->DistributeDofs(analysis_type);
    29 
     28                nodes->DistributeDofs(analysis_type,GsetEnum);
     29                nodes->DistributeDofs(analysis_type,FsetEnum);
     30                nodes->DistributeDofs(analysis_type,SsetEnum);
    3031        }
    3132

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

    r4217 r5772  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void SpcNodesx(Vec* pyg, Nodes* nodes,Constraints* constraints,int analysis_type){
     12void SpcNodesx(Nodes* nodes,Constraints* constraints,int analysis_type){
    1313
    1414        int i;
    15         int numberofdofs;
    16         int gsize;
     15        Node* node=NULL;
     16        int nodeid;
     17        int dof;
     18        double value;
    1719
    18         /*output: */
    19         Vec yg=NULL;
     20        for(i=0;i<constraints->Size();i++){
    2021
    21         /*First, recover number of dofs from nodes: */
    22         numberofdofs=nodes->NumberOfDofs(analysis_type);
     22                Object* object=(Object*)constraints->GetObjectByOffset(i);
    2323
    24         if(numberofdofs){
    25                
    26                 /*Allocate yg: */
    27                 yg=NewVec(numberofdofs);
     24                /*Check constraints is a single point constraint (spc): */
     25                if(object->Enum()==SpcEnum){
    2826
    29                 /*Now, go through constraints, and update the nodes and the constraint vector at the same time: */
    30                 constraints->SetupSpcs(nodes,yg,analysis_type);
     27                        Spc* spc=(Spc*)object;
    3128
    32                 /*Specify numentries: */
    33                 VecGetSize(yg,&gsize);
     29                        if(spc->InAnalysis(analysis_type)){
     30
     31                                /*Ok, constraints object is a constraint. Get the nodeid from the node it applies to: */
     32                                nodeid=spc->GetNodeId();
     33                                dof=spc->GetDof();
     34                                value=spc->GetValue();
     35
     36                                /*Now, chase through nodes and find the corect node: */
     37                                node=(Node*)nodes->GetObjectById(NULL,nodeid);
     38
     39                                /*Apply constraint: */
     40                                if(node){ //in case the spc is dealing with a node on another cpu
     41                                        node->ApplyConstraint(dof,value);
     42                                }
     43                        }
     44
     45                }
    3446        }
    3547
    36         /*Assign output pointers: */
    37         *pyg=yg;
    3848}

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

    r4236 r5772  
    1111
    1212/* local prototypes: */
    13 void SpcNodesx(Vec* pyg, Nodes* nodesin,Constraints* constraints,int analysis_type);
     13void SpcNodesx(Nodes* nodesin,Constraints* constraints,int analysis_type);
    1414
    1515#endif  /* _SPCNODESX_H */

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

    r5703 r5772  
    99#include "../../EnumDefinitions/EnumDefinitions.h"
    1010
    11 void 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){
     11void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){
    1212       
    1313        /*intermediary: */
    14         int      i,j,gsize;
     14        int      i,j;
     15        int      gsize,fsize,ssize;
    1516        int      connectivity, numberofdofspernode;
    1617        int      analysis_type, configuration_type;
    1718        Element *element = NULL;
    1819        Load    *load    = NULL;
     20        bool    buildkff=false;
    1921       
    2022        /*output: */
    2123        Mat    Kgg  = NULL;
     24        Mat    Kff  = NULL;
     25        Mat    Kfs  = NULL;
    2226        Vec    pg   = NULL;
     27        Vec    pf   = NULL;
    2328        double kmax = 0;
    2429
    … …  
    3136        parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    3237        parameters->FindParam(&connectivity,ConnectivityEnum);
     38        parameters->FindParam(&buildkff,KffEnum);
    3339
    3440        /*Get size of matrix: */
    35         gsize=nodes->NumberOfDofs(configuration_type);
    36         numberofdofspernode=nodes->MaxNumDofs(configuration_type);
     41        gsize=nodes->NumberOfDofs(configuration_type,GsetEnum);
     42        fsize=nodes->NumberOfDofs(configuration_type,FsetEnum);
     43        ssize=nodes->NumberOfDofs(configuration_type,SsetEnum);
     44        numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum);
    3745
    3846        /*Checks in debugging mode {{{1*/
    … …  
    4654        if(kflag){
    4755
    48                 Kgg=NewMat(gsize,gsize,NULL,&connectivity,&numberofdofspernode);
     56                if(!buildkff)Kgg=NewMat(gsize,gsize,NULL,&connectivity,&numberofdofspernode);
     57                else{
     58                        Kff=NewMat(fsize,fsize,NULL,&connectivity,&numberofdofspernode);
     59                        Kfs=NewMat(fsize,ssize,NULL,&connectivity,&numberofdofspernode);
     60                }
    4961
    5062                /*Fill stiffness matrix from elements: */
    5163                for (i=0;i<elements->Size();i++){
    5264                        element=(Element*)elements->GetObjectByOffset(i);
    53                         element->CreateKMatrix(Kgg);
     65                        element->CreateKMatrix(Kgg,Kff,Kfs);
    5466                }
    5567
    … …  
    5769                for (i=0;i<loads->Size();i++){
    5870                        load=(Load*)loads->GetObjectByOffset(i);
    59                         if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kgg);
     71                        if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kgg,Kff,Kfs);
    6072                }
    6173
    6274                /*Assemble matrix and compress matrix to save memory: */
    63                 MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY);
    64                 MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY);
    65                 MatCompress(Kgg);
     75                if(!buildkff){
     76                        MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY);
     77                        MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY);
     78                        MatCompress(Kgg);
     79                }
     80                else{
     81                        MatAssemblyBegin(Kff,MAT_FINAL_ASSEMBLY);
     82                        MatAssemblyEnd(Kff,MAT_FINAL_ASSEMBLY);
     83                        MatCompress(Kff);
     84
     85                        MatAssemblyBegin(Kfs,MAT_FINAL_ASSEMBLY);
     86                        MatAssemblyEnd(Kfs,MAT_FINAL_ASSEMBLY);
     87                        MatCompress(Kfs);
     88                }
    6689        }
    6790        if(pflag){
    … …  
    7295                for (i=0;i<elements->Size();i++){
    7396                        element=(Element*)elements->GetObjectByOffset(i);
    74                         element->CreatePVector(pg);
     97                        element->CreatePVector(pg,pf);
    7598                }
    7699
    … …  
    78101                for (i=0;i<loads->Size();i++){
    79102                        load=(Load*)loads->GetObjectByOffset(i);
    80                         if (load->InAnalysis(configuration_type)) load->CreatePVector(pg);
     103                        if (load->InAnalysis(configuration_type)) load->CreatePVector(pg,pf);
    81104                }
    82105
    83                 VecAssemblyBegin(pg);
    84                 VecAssemblyEnd(pg);
     106                if(!buildkff){
     107                        VecAssemblyBegin(pg);
     108                        VecAssemblyEnd(pg);
     109                }
     110                else{
     111                        VecAssemblyBegin(pf);
     112                        VecAssemblyEnd(pf);
     113                }
    85114        }
    86115
     116       
    87117        /*Now, deal with penalties*/
    88118        if(penalty_kflag){
    89119
    90120                /*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */
    91                 MatNorm(Kgg,NORM_INFINITY,&kmax);
     121                if(!buildkff)MatNorm(Kgg,NORM_INFINITY,&kmax);
     122                else MatNorm(Kff,NORM_INFINITY,&kmax);
    92123
    93124                /*Fill stiffness matrix from loads: */
    94125                for (i=0;i<loads->Size();i++){
    95126                        load=(Load*)loads->GetObjectByOffset(i);
    96                         if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kgg,kmax);
     127                        if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kgg,Kff,Kfs,kmax);
    97128                }
    98129
    99130                /*Assemble matrix and compress matrix to save memory: */
    100                 MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY);
    101                 MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY);
    102                 MatCompress(Kgg);
     131                if(!buildkff){
     132                        MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY);
     133                        MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY);
     134                        MatCompress(Kgg);
     135                }
     136                else{
     137                        MatAssemblyBegin(Kff,MAT_FINAL_ASSEMBLY);
     138                        MatAssemblyEnd(Kff,MAT_FINAL_ASSEMBLY);
     139                        MatCompress(Kff);
     140
     141                        MatAssemblyBegin(Kfs,MAT_FINAL_ASSEMBLY);
     142                        MatAssemblyEnd(Kfs,MAT_FINAL_ASSEMBLY);
     143                        MatCompress(Kfs);
     144                }
    103145        }
    104146        if(penalty_pflag){
    … …  
    107149                for (i=0;i<loads->Size();i++){
    108150                        load=(Load*)loads->GetObjectByOffset(i);
    109                         if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pg,kmax);
     151                        if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pg,pf,kmax);
    110152                }
    111153
    112                 VecAssemblyBegin(pg);
    113                 VecAssemblyEnd(pg);
     154                if(!buildkff){
     155                        VecAssemblyBegin(pg);
     156                        VecAssemblyEnd(pg);
     157                }
     158                else{
     159                        VecAssemblyBegin(pf);
     160                        VecAssemblyEnd(pf);
     161                }
    114162        }
    115163       
    116164        /*Assign output pointers: */
    117         *pKgg=Kgg;
    118         *ppg=pg;
     165        if(pKgg) *pKgg=Kgg;
     166        if(pKff) *pKff=Kff;
     167        if(pKfs) *pKfs=Kfs;
     168        if(ppg)  *ppg=pg;
     169        if(ppf)  *ppf=pf;
    119170        if(pkmax) *pkmax=kmax;
    120171}

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

    r5693 r5772  
    1010
    1111/* local prototypes: */
    12 void SystemMatricesx(Mat* pKgg, Vec* ppg,double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,
     12void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,
    1313                        bool kflag=true,bool pflag=true,bool penalty_kflag=true,bool penalty_pflag=true);
    1414

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

    r5709 r5772  
    2020#include "./ContourToNodesx/ContourToNodesx.h"
    2121#include "./CostFunctionx/CostFunctionx.h"
     22#include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h"
    2223#include "./DakotaResponsesx/DakotaResponsesx.h"
    2324#include "./ElementConnectivityx/ElementConnectivityx.h"
    … …  
    7273#include "./Qmux/Qmux.h"
    7374#include "./Reduceloadfromgtofx/Reduceloadfromgtofx.h"
     75#include "./Reduceloadx/Reduceloadx.h"
    7476#include "./Reducematrixfromgtofx/Reducematrixfromgtofx.h"
    7577#include "./Reducevectorgtosx/Reducevectorgtosx.h"

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

    r5254 r5772  
    2121DofIndexing::DofIndexing(){
    2222
    23         int i;
    24         this->numberofdofs=UNDEF;
     23        this->gsize=UNDEF;
     24        this->fsize=UNDEF;
     25        this->ssize=UNDEF;
    2526        this->clone=0;
    2627        this->f_set=NULL;
    2728        this->s_set=NULL;
    28         this->doflist=NULL;
     29        this->svalues=NULL;
    2930        this->doftype=NULL;
    30 }
    31 /*}}}*/
    32 /*FUNCTION DofIndexing::DofIndexing(int numberofdofs){{{1*/
    33 DofIndexing::DofIndexing(int in_numberofdofs){
    34         this->Init(in_numberofdofs,NULL);
     31        this->gdoflist=NULL;
     32        this->fdoflist=NULL;
     33        this->sdoflist=NULL;
     34
     35}
     36/*}}}*/
     37/*FUNCTION DofIndexing::DofIndexing(int gsize){{{1*/
     38DofIndexing::DofIndexing(int in_gsize){
     39        this->Init(in_gsize,NULL);
    3540}
    3641/*}}}*/
    … …  
    3944
    4045        int i;
    41         this->numberofdofs=in->numberofdofs;
     46        this->gsize=in->gsize;
     47        this->fsize=in->fsize;
     48        this->ssize=in->ssize;
     49       
    4250        this->clone=in->clone;
    4351
    44         this->f_set=(int*)xmalloc(this->numberofdofs*sizeof(int));
    45         this->s_set=(int*)xmalloc(this->numberofdofs*sizeof(int));
    46         this->doflist=(int*)xmalloc(this->numberofdofs*sizeof(int));
    47         if (in->doftype) this->doftype=(int*)xmalloc(this->numberofdofs*sizeof(int));
    48         else this->doftype=NULL;
    49 
    50         for(i=0;i<this->numberofdofs;i++){
    51                 this->f_set[i]=in->f_set[i];
    52                 this->s_set[i]=in->s_set[i];
    53                 this->doflist[i]=in->doflist[i];
    54                 if (in->doftype) this->doftype[i]=in->doftype[i];
    55         }
     52        if(this->gsize){
     53                this->f_set=(bool*)xmalloc(this->gsize*sizeof(bool));
     54                this->s_set=(bool*)xmalloc(this->gsize*sizeof(bool));
     55                this->svalues=(double*)xmalloc(this->gsize*sizeof(int));
     56                if(in->doftype)this->doftype=(int*)xmalloc(this->gsize*sizeof(int));
     57                this->gdoflist=(int*)xmalloc(this->gsize*sizeof(int));
     58        }
     59        else{
     60                this->f_set=NULL;
     61                this->s_set=NULL;
     62                this->svalues=NULL;
     63                this->doftype=NULL;
     64                this->gdoflist=NULL;
     65        }
     66        if(this->fsize)this->fdoflist=(int*)xmalloc(this->fsize*sizeof(int)); else this->fdoflist=NULL;
     67        if(this->ssize)this->sdoflist=(int*)xmalloc(this->ssize*sizeof(int)); else this->sdoflist=NULL;
     68
     69        if(this->gsize){
     70                memcpy(this->f_set,in->f_set,this->gsize*sizeof(bool));
     71                memcpy(this->s_set,in->s_set,this->gsize*sizeof(bool));
     72                memcpy(this->svalues,in->svalues,this->gsize*sizeof(double));
     73                if(this->doftype)memcpy(this->doftype,in->doftype,this->gsize*sizeof(int));
     74                memcpy(this->gdoflist,in->gdoflist,this->gsize*sizeof(int));
     75        }
     76        if(this->fsize)memcpy(this->fdoflist,in->fdoflist,this->fsize*sizeof(int));
     77        if(this->ssize)memcpy(this->sdoflist,in->sdoflist,this->ssize*sizeof(int));
     78
    5679}
    5780/*}}}*/
    … …  
    6184        xfree((void**)&f_set);
    6285        xfree((void**)&s_set);
    63         xfree((void**)&doflist);
     86        xfree((void**)&svalues);
    6487        xfree((void**)&doftype);
     88        xfree((void**)&gdoflist);
     89        xfree((void**)&fdoflist);
     90        xfree((void**)&sdoflist);
    6591
    6692}
    6793/*}}}*/
    6894/*FUNCTION DofIndexing::Init{{{1*/
    69 void DofIndexing::Init(int in_numberofdofs,int* in_doftype){
    70 
    71         int i;
    72         this->numberofdofs=in_numberofdofs;
     95void DofIndexing::Init(int in_gsize,int* in_doftype){
     96
     97        int i;
     98        this->gsize=in_gsize;
     99       
    73100        this->clone=0;
    74101
    75102        /*allocate: */
    76         this->f_set=(int*)xmalloc(this->numberofdofs*sizeof(int));
    77         this->s_set=(int*)xmalloc(this->numberofdofs*sizeof(int));
    78         this->doflist=(int*)xmalloc(this->numberofdofs*sizeof(int));
    79         if (in_doftype) this->doftype=(int*)xmalloc(this->numberofdofs*sizeof(int));
    80 
    81         for (i=0;i<this->numberofdofs;i++){
     103        if(this->gsize){
     104                this->f_set=(bool*)xmalloc(this->gsize*sizeof(bool));
     105                this->s_set=(bool*)xmalloc(this->gsize*sizeof(bool));
     106                this->svalues=(double*)xmalloc(this->gsize*sizeof(double));
     107                if(in_doftype)this->doftype=(int*)xmalloc(this->gsize*sizeof(int));
     108                this->gdoflist=(int*)xmalloc(this->gsize*sizeof(int));
     109        }
     110
     111        for (i=0;i<this->gsize;i++){
    82112                /*assume dof is free, no constraints, no rigid body constraint: */
    83                 this->f_set[i]=1;
    84                 this->s_set[i]=0;
    85                 this->doflist[i]=UNDEF;
    86                 if(in_doftype) this->doftype[i]=in_doftype[i];
    87         }
     113                this->f_set[i]=true;
     114                this->s_set[i]=false;
     115                if(this->doftype)this->doftype[i]=in_doftype[i];
     116                this->svalues[i]=0; //0 constraint is the default value
     117        }
     118}
     119/*}}}*/
     120/*FUNCTION DofIndexing::InitSet{{{1*/
     121void DofIndexing::InitSet(int setenum){
     122
     123        int i;
     124        int size=0;
     125
     126        /*go through sets, and figure out how many dofs belong to this set, except for g-set,
     127         * which has already been initialized: */
     128        if(setenum==FsetEnum){
     129                size=0;
     130                for(i=0;i<this->gsize;i++) if(f_set[i])size++;
     131                this->fsize=size;
     132                if(this->fsize)this->fdoflist=(int*)xmalloc(size*sizeof(int));
     133                else this->fdoflist=NULL;
     134        }
     135        else if(setenum==SsetEnum){
     136                size=0;
     137                for(i=0;i<this->gsize;i++) if(s_set[i])size++;
     138                this->ssize=size;
     139                if(this->ssize)this->sdoflist=(int*)xmalloc(size*sizeof(int));
     140                else this->sdoflist=NULL;
     141        }
     142        else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
     143
     144
    88145}
    89146/*}}}*/
    … …  
    96153
    97154        printf("DofIndexing:\n");
    98         printf("   numberofdofs: %i\n",numberofdofs);
     155        printf("   gsize: %i\n",gsize);
    99156        printf("   clone: %i\n",clone);
    100157}
    … …  
    106163
    107164        printf("DofIndexing:\n");
    108         printf("   numberofdofs: %i\n",numberofdofs);
     165        printf("   gsize: %i\n",gsize);
     166        printf("   fsize: %i\n",fsize);
     167        printf("   ssize: %i\n",ssize);
    109168        printf("   clone: %i\n",clone);
    110169       
    111170        printf("   set membership: f,s sets \n");
    112         for(i=0;i<numberofdofs;i++){
    113                 printf("      dof %i: %i %i\n",i,f_set[i],s_set[i]);
    114         }
    115         printf("\n");
    116 
    117         printf("   doflist: |");
    118         for(i=0;i<numberofdofs;i++){
    119                 printf(" %i |",doflist[i]);
     171        for(i=0;i<gsize;i++){
     172                printf("      dof %i: %s %s\n",i,f_set[i]?"true":"false",s_set[i]?"true":"false");
     173        }
     174
     175        printf("   svalues (%i): |",this->ssize);
     176        for(i=0;i<this->gsize;i++){
     177                if(this->s_set[i])printf(" %g |",svalues[i]);
    120178        }
    121179        printf("\n");
    … …  
    123181        if(doftype){
    124182                printf("   doftype: |");
    125                 for(i=0;i<numberofdofs;i++){
     183                for(i=0;i<gsize;i++){
    126184                        printf(" %i |",doftype[i]);
    127185                }
    128186                printf("\n");
    129187        }
    130         else printf("   NULL doftype");
    131 
     188        else printf("   doftype: NULL\n");
     189
     190        printf("   g_doflist (%i): |",this->gsize);
     191        for(i=0;i<this->gsize;i++){
     192                printf(" %i |",gdoflist[i]);
     193        }
     194        printf("\n");
     195
     196        printf("   f_doflist (%i): |",this->fsize);
     197        for(i=0;i<this->fsize;i++){
     198                printf(" %i |",fdoflist[i]);
     199        }
     200        printf("\n");
     201
     202        printf("   s_doflist (%i): |",this->ssize);
     203        for(i=0;i<this->ssize;i++){
     204                printf(" %i |",sdoflist[i]);
     205        }
     206        printf("\n");
    132207}               
    133208/*}}}*/
    … …  
    145220        memcpy(&enum_type,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int);
    146221
    147         memcpy(&numberofdofs,marshalled_dataset,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs);
     222        /*easy part: */
     223        memcpy(&gsize,marshalled_dataset,sizeof(gsize));marshalled_dataset+=sizeof(gsize);
     224        memcpy(&fsize,marshalled_dataset,sizeof(fsize));marshalled_dataset+=sizeof(fsize);
     225        memcpy(&ssize,marshalled_dataset,sizeof(ssize));marshalled_dataset+=sizeof(ssize);
     226        memcpy(&flagdoftype,marshalled_dataset,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype);
    148227        memcpy(&clone,marshalled_dataset,sizeof(clone));marshalled_dataset+=sizeof(clone);
    149228       
    150229        /*Allocate: */
    151         this->f_set=(int*)xmalloc(this->numberofdofs*sizeof(int));
    152         this->s_set=(int*)xmalloc(this->numberofdofs*sizeof(int));
    153         this->doflist=(int*)xmalloc(this->numberofdofs*sizeof(int));
    154 
    155         memcpy(f_set,marshalled_dataset,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int);
    156         memcpy(s_set,marshalled_dataset,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int);
    157         memcpy(doflist,marshalled_dataset,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int);
    158        
    159         memcpy(&flagdoftype,marshalled_dataset,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype);
    160         if(flagdoftype){ // there is a hook so demarshall it
    161                 this->doftype=(int*)xmalloc(this->numberofdofs*sizeof(int));
    162                 memcpy(doftype,marshalled_dataset,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int);
    163         }
    164         else this->doftype=NULL; //There is no coupling, so doftype is NULL
     230        if(this->gsize){
     231                this->f_set=(bool*)xmalloc(this->gsize*sizeof(bool));
     232                this->s_set=(bool*)xmalloc(this->gsize*sizeof(bool));
     233                this->svalues=(double*)xmalloc(this->ssize*sizeof(double));
     234                if(flagdoftype)this->doftype=(int*)xmalloc(this->gsize*sizeof(int));
     235                else           this->doftype=NULL;
     236                this->gdoflist=(int*)xmalloc(this->gsize*sizeof(int));
     237        }
     238        else{
     239                this->f_set=NULL;
     240                this->s_set=NULL;
     241                this->svalues=NULL;
     242                this->doftype=NULL;
     243                this->gdoflist=NULL;
     244        }
     245        if(this->fsize)this->fdoflist=(int*)xmalloc(this->fsize*sizeof(int));
     246        else           this->fdoflist=NULL;
     247        if(this->ssize)this->sdoflist=(int*)xmalloc(this->ssize*sizeof(int));
     248        else           this->sdoflist=NULL;
     249
     250        /*Copy arrays: */
     251        if(this->gsize){
     252                memcpy(f_set,marshalled_dataset,gsize*sizeof(bool));marshalled_dataset+=gsize*sizeof(bool);
     253                memcpy(s_set,marshalled_dataset,gsize*sizeof(bool));marshalled_dataset+=gsize*sizeof(bool);
     254                memcpy(this->svalues,marshalled_dataset,this->gsize*sizeof(double));marshalled_dataset+=this->gsize*sizeof(double);
     255                if(flagdoftype){ memcpy(doftype,marshalled_dataset,gsize*sizeof(int));marshalled_dataset+=gsize*sizeof(int); }
     256                memcpy(this->gdoflist,marshalled_dataset,this->gsize*sizeof(int));marshalled_dataset+=this->gsize*sizeof(int);
     257        }
     258       
     259        if(this->fsize){ memcpy(this->fdoflist,marshalled_dataset,this->fsize*sizeof(int));marshalled_dataset+=this->fsize*sizeof(int); }
     260        if(this->ssize){ memcpy(this->sdoflist,marshalled_dataset,this->ssize*sizeof(int));marshalled_dataset+=this->ssize*sizeof(int); }
    165261
    166262        /*return: */
    … …  
    174270        char* marshalled_dataset=NULL;
    175271        int   enum_type=0;
    176         int   flagdoftype; //to indicate if there are some doftype or if NULL
     272        bool  flagdoftype; //to indicate if there are some doftype or if NULL
    177273
    178274        /*recover marshalled_dataset: */
    179275        marshalled_dataset=*pmarshalled_dataset;
    180276
     277        /*preliminary: */
     278        if(this->doftype)flagdoftype=true;
     279        else             flagdoftype=false;
     280
    181281        /*get enum type of DofIndexing: */
    182282        enum_type=DofIndexingEnum;
    … …  
    186286       
    187287        /*marshall DofIndexing data: */
    188         memcpy(marshalled_dataset,&numberofdofs,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs);
     288        memcpy(marshalled_dataset,&gsize,sizeof(gsize));marshalled_dataset+=sizeof(gsize);
     289        memcpy(marshalled_dataset,&fsize,sizeof(fsize));marshalled_dataset+=sizeof(fsize);
     290        memcpy(marshalled_dataset,&ssize,sizeof(ssize));marshalled_dataset+=sizeof(ssize);
     291        memcpy(marshalled_dataset,&flagdoftype,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype);
    189292        memcpy(marshalled_dataset,&clone,sizeof(clone));marshalled_dataset+=sizeof(clone);
    190         memcpy(marshalled_dataset,f_set,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int);
    191         memcpy(marshalled_dataset,s_set,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int);
    192         memcpy(marshalled_dataset,doflist,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int);
    193 
    194         /*marshall doftype if there are some and add a flag*/
    195         if(this->doftype){
    196                 flagdoftype=1;
    197                 memcpy(marshalled_dataset,&flagdoftype,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype);
    198                 memcpy(marshalled_dataset,doftype,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int);
    199         }
    200         else{
    201                 flagdoftype=0;
    202                 memcpy(marshalled_dataset,&flagdoftype,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype);
    203         }
     293       
     294       
     295        if(this->gsize){
     296                memcpy(marshalled_dataset,f_set,gsize*sizeof(bool));marshalled_dataset+=gsize*sizeof(bool);
     297                memcpy(marshalled_dataset,s_set,gsize*sizeof(bool));marshalled_dataset+=gsize*sizeof(bool);
     298                memcpy(marshalled_dataset,svalues,gsize*sizeof(double)); marshalled_dataset+=gsize*sizeof(double);
     299                if(flagdoftype){ memcpy(marshalled_dataset,doftype,gsize*sizeof(int)); marshalled_dataset+=gsize*sizeof(int); }
     300                memcpy(marshalled_dataset,gdoflist,gsize*sizeof(int)); marshalled_dataset+=gsize*sizeof(int);
     301        }
     302        if(this->fsize){ memcpy(marshalled_dataset,fdoflist,fsize*sizeof(int)); marshalled_dataset+=fsize*sizeof(int);}
     303        if(this->ssize){ memcpy(marshalled_dataset,sdoflist,ssize*sizeof(int)); marshalled_dataset+=ssize*sizeof(int);}
    204304
    205305        *pmarshalled_dataset=marshalled_dataset;
    … …  
    210310int   DofIndexing::MarshallSize(){
    211311
    212         int doftype_size=0;
    213        
    214         /*FInd size of doftype*/
    215         doftype_size+=sizeof(int); //Flag 0 or 1
    216         if (this->doftype){
    217                 doftype_size+=numberofdofs*sizeof(int);
    218         }
    219 
    220         return sizeof(numberofdofs)+
    221                 sizeof(clone)+
    222                 numberofdofs*sizeof(int)+
    223                 numberofdofs*sizeof(int)+
    224                 numberofdofs*sizeof(int)+
    225                 doftype_size+
    226                 sizeof(int); //sizeof(int) for enum type
    227 }
    228 /*}}}*/
     312        int size=0;
     313
     314        size+=4*sizeof(int)+sizeof(bool);
     315        if(this->gsize){
     316                size+= 2*this->gsize*sizeof(bool)+
     317                           this->gsize*sizeof(double)+
     318                           this->gsize*sizeof(int);
     319                if(this->doftype)size+=this->gsize*sizeof(int);
     320        }
     321        if(this->fsize)size+=this->fsize*sizeof(int);
     322        if(this->ssize)size+=this->ssize*sizeof(int);
     323
     324        size+=sizeof(int); //sizeof(int) for enum type
     325
     326        return size;
     327}
     328/*}}}*/

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

    r5254 r5772  
    1010        public:
    1111
    12                 int numberofdofs; //number of dofs for a node, variable
     12                /*sizes: */
     13                int gsize; //number of dofs for a node
     14                int fsize; //number of dofs solver for
     15                int ssize; //number of constrained dofs
    1316
    1417                /*partitioning: */
    … …  
    1619
    1720                /*boundary conditions sets: */
    18                 int*     f_set; //set on which we solve
    19                 int*     s_set; //set on which boundary conditions (dirichlet) are applied
     21                bool*     f_set; //is dof on f-set (on which we solve)
     22                bool*     s_set; //is dof on s-set (on which boundary conditions -dirichlet- are applied)
     23                double*   svalues; //list of constraint values. size g_size, for ease of use.
    2024
     25                /*types of dofs: */
     26                int*     doftype; //approximation type of the dofs (used only for coupling), size g_size
     27               
    2128                /*list of degrees of freedom: */
    22                 int*     doflist; //dof list on which we solve
    23                 int*     doftype; //approximation type of the dofs (used only for coupling)
     29                int*     gdoflist; //dof list in g_set
     30                int*     fdoflist; //dof list in f_set
     31                int*     sdoflist; //dof list in s_set
     32
    2433
    2534                /*DofIndexing constructors, destructors {{{1*/
    2635                DofIndexing();
    27                 DofIndexing(int numberofdofs);
    28                 void Init(int numberofdofs,int* doftype);
     36                DofIndexing(int g_size);
     37                void Init(int g_size,int* doftype);
     38                void InitSet(int setenum);
    2939                DofIndexing(DofIndexing* properties);
    3040                ~DofIndexing();

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

    r5745 r5772  
    2828                virtual void   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0;
    2929                virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0;
    30                 virtual void   CreateKMatrix(Mat Kgg)=0;
    31                 virtual void   CreatePVector(Vec pg)=0;
     30                virtual void   CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs)=0;
     31                virtual void   CreatePVector(Vec pg, Vec pf)=0;
    3232                virtual void   GetSolutionFromInputs(Vec solution)=0;
    3333                virtual int    GetNodeIndex(Node* node)=0;

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

    r5769 r5772  
    367367        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    368368        if (analysis_type==DiagnosticHorizAnalysisEnum){
    369                 InputUpdateFromSolutionDiagnosticHoriz( solution);
     369                int approximation;
     370                inputs->GetParameterValue(&approximation,ApproximationEnum);
     371                if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
     372                        InputUpdateFromSolutionDiagnosticStokes( solution);
     373                }
     374                else{
     375                        InputUpdateFromSolutionDiagnosticHoriz( solution);
     376                }
    370377        }
    371378        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    685692/*}}}*/
    686693/*FUNCTION Penta::CreateKMatrix {{{1*/
    687 void  Penta::CreateKMatrix(Mat Kgg){
     694void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
    688695
    689696        int analysis_type;
    … …  
    773780/*}}}*/
    774781/*FUNCTION Penta::CreatePVector {{{1*/
    775 void  Penta::CreatePVector(Vec pg){
     782void  Penta::CreatePVector(Vec pg, Vec pf){
    776783
    777784        int analysis_type;
    … …  
    894901                        GetSolutionFromInputsDiagnosticStokes(solution);
    895902                }
    896                 else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum ||approximation==HutterApproximationEnum){
     903                else{
    897904                        GetSolutionFromInputsDiagnosticHoriz(solution);
    898905                }
    899                 else if(approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum){
    900                         return; //do nothing the elements around with update the solution
    901                 }
    902                 else ISSMERROR("approximation type : %i (%s) not implemented yet",approximation,EnumToString(approximation));
    903906        }
    904907        else if(analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    18921895                        this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum));
    18931896                }
    1894                 else if (*(iomodel->elements_type+index)==PattynStokesApproximationEnum){
    1895                         this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum));
    1896                 }
    18971897                else if (*(iomodel->elements_type+index)==NoneApproximationEnum){
    18981898                        this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum));
    … …  
    19881988        /*Spawn Tria element from the base of the Penta: */
    19891989        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1990         tria->CreateKMatrix(Kgg);
     1990        tria->CreateKMatrix(Kgg,NULL,NULL);
    19911991        delete tria->matice; delete tria;
    19921992
    … …  
    20182018        /*Spawn Tria element from the base of the Penta: */
    20192019        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2020         tria->CreateKMatrix(Kgg);
     2020        tria->CreateKMatrix(Kgg,NULL,NULL);
    20212021        delete tria->matice; delete tria;
    20222022
    … …  
    20962096        /* Get node coordinates and dof list: */
    20972097        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESp);
    2098         tria->GetDofList(&doflistm,MacAyealApproximationEnum);  //Pattyn dof list
    2099         GetDofList(&doflistp,PattynApproximationEnum); //MacAyeal dof list
     2098        tria->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);  //Pattyn dof list
     2099        GetDofList(&doflistp,PattynApproximationEnum,GsetEnum); //MacAyeal dof list
    21002100
    21012101        /*Retrieve all inputs we will be needing: */
    … …  
    21832183        if(IsOnWater())return;
    21842184
    2185         GetDofList(&doflist);
     2185        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    21862186
    21872187        /*Spawn 3 beam elements: */
    … …  
    23212321                        /*Call Tria function*/
    23222322                        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2323                         tria->CreateKMatrix(Kgg);
     2323                        tria->CreateKMatrix(Kgg,NULL,NULL);
    23242324                        delete tria->matice; delete tria;
    23252325
    … …  
    23432343                /* Get node coordinates and dof list: */
    23442344                GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES);
    2345                 tria->GetDofList(&doflist,MacAyealApproximationEnum);  //Pattyn dof list
     2345                tria->GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum);  //Pattyn dof list
    23462346
    23472347                /*Retrieve all inputs we will be needing: */
    … …  
    23992399                         * grids: */
    24002400
    2401                         tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg);
     2401                        tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg,NULL,NULL);
    24022402                }
    24032403
    … …  
    24762476        /* Get node coordinates and dof list: */
    24772477        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2478         GetDofList(&doflist,PattynApproximationEnum);
     2478        GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    24792479
    24802480        /*Retrieve all inputs we will be needing: */
    … …  
    25962596        /* Get node coordinates and dof list: */
    25972597        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2598         GetDofList(&doflist,StokesApproximationEnum);
     2598        GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
    25992599
    26002600        /*Retrieve all inputs we will be needing: */
    … …  
    27562756        /* Get node coordinates and dof list: */
    27572757        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2758         GetDofList(&doflist);
     2758        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    27592759
    27602760        /* Start  looping on the number of gaussian points: */
    … …  
    28272827        /*Spawn Tria element from the base of the Penta: */
    28282828        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2829         tria->CreateKMatrix(Kgg);
     2829        tria->CreateKMatrix(Kgg,NULL,NULL);
    28302830        delete tria->matice; delete tria;
    28312831
    … …  
    28522852        /*Spawn Tria element from the base of the Penta: */
    28532853        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2854         tria->CreateKMatrix(Kgg);
     2854        tria->CreateKMatrix(Kgg,NULL,NULL);
    28552855        delete tria->matice; delete tria;
    28562856        return;
    … …  
    29202920        /* Get node coordinates and dof list: */
    29212921        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2922         GetDofList(&doflist);
     2922        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    29232923
    29242924        // /*recovre material parameters: */
    … …  
    31093109        /*Spawn Tria element from the base of the Penta: */
    31103110        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3111         tria->CreatePVector(pg);
     3111        tria->CreatePVector(pg,NULL);
    31123112        delete tria->matice; delete tria;
    31133113
    … …  
    31373137        /*Spawn Tria element from the base of the Penta: */
    31383138        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3139         tria->CreatePVector(pg);
     3139        tria->CreatePVector(pg,NULL);
    31403140        delete tria->matice; delete tria;
    31413141
    … …  
    31783178
    31793179        /*recover doflist: */
    3180         GetDofList(&doflist);
     3180        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    31813181
    31823182        /*recover parameters: */
    … …  
    32753275                 *and use its CreatePVector functionality to return an elementary load vector: */
    32763276                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3277                 tria->CreatePVector(pg);
     3277                tria->CreatePVector(pg,NULL);
    32783278                delete tria->matice; delete tria;
    32793279                return;
    … …  
    33173317        /* Get node coordinates and dof list: */
    33183318        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3319         GetDofList(&doflist,PattynApproximationEnum);
     3319        GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    33203320
    33213321        /*Retrieve all inputs we will be needing: */
    … …  
    34293429        /* Get node coordinates and dof list: */
    34303430        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3431         GetDofList(&doflist,StokesApproximationEnum);
     3431        GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
    34323432
    34333433        /*Retrieve all inputs we will be needing: */
    … …  
    35933593        /* Get node coordinates and dof list: */
    35943594        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3595         GetDofList(&doflist);
     3595        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    35963596
    35973597        /*Retrieve all inputs we will be needing: */
    … …  
    36563656        /*Spawn Tria element from the base of the Penta: */
    36573657        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3658         tria->CreatePVector(pg);
     3658        tria->CreatePVector(pg,NULL);
    36593659        delete tria->matice; delete tria;
    36603660
    … …  
    36803680        /*Spawn Tria element from the base of the Penta: */
    36813681        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3682         tria->CreatePVector(pg);
     3682        tria->CreatePVector(pg,NULL);
    36833683        delete tria->matice; delete tria;
    36843684        return;
    … …  
    37443744        /* Get node coordinates and dof list: */
    37453745        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3746         GetDofList(&doflist);
     3746        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    37473747
    37483748        /*recovre material parameters: */
    … …  
    38213821/*}}}*/
    38223822/*FUNCTION Penta::GetDofList {{{1*/
    3823 void  Penta::GetDofList(int** pdoflist,int approximation_enum){
     3823void  Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){
    38243824
    38253825        int i,j;
    … …  
    38323832        /*First, figure out size of doflist: */
    38333833        for(i=0;i<6;i++){
    3834                 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum);
     3834                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    38353835        }
    38363836
    … …  
    38413841        count=0;
    38423842        for(i=0;i<6;i++){
    3843                 nodes[i]->GetDofList(doflist+count,approximation_enum);
    3844                 count+=nodes[i]->GetNumberOfDofs(approximation_enum);
     3843                nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
     3844                count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    38453845        }
    38463846
    … …  
    39953995        /*If the element is a coupling, do nothing: every grid is also on an other elements
    39963996         * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
    3997         GetDofList(&doflist,approximation);
     3997        if(approximation==MacAyealPattynApproximationEnum){
     3998         return;
     3999        }
     4000        else{
     4001                GetDofList(&doflist,approximation,GsetEnum);
     4002        }
    39984003
    39994004        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    40324037
    40334038        /*Get dof list: */
    4034         GetDofList(&doflist);
     4039        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    40354040        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
    40364041        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    … …  
    40704075
    40714076        /*Get dof list: */
    4072         GetDofList(&doflist);
     4077        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    40734078        Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
    40744079
    … …  
    41064111
    41074112        /*Get dof list: */
    4108         GetDofList(&doflist);
     4113        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    41094114        Input* vx_input=inputs->GetInput(VxEnum);       ISSMASSERT(vx_input);
    41104115        Input* vy_input=inputs->GetInput(VyEnum);       ISSMASSERT(vy_input);
    … …  
    41524157
    41534158        /*Get dof list: */
    4154         GetDofList(&doflist);
     4159        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    41554160        Input* t_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(t_input);
    41564161
    … …  
    44424447                InputUpdateFromSolutionDiagnosticPattyn(solution);
    44434448        }
    4444         else if (approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
    4445                 InputUpdateFromSolutionDiagnosticStokes(solution);
    4446         }
    44474449        else if (approximation==MacAyealPattynApproximationEnum){
    44484450                InputUpdateFromSolutionDiagnosticMacAyealPattyn(solution);
    4449         }
    4450         else if (approximation==PattynStokesApproximationEnum){
    4451                 InputUpdateFromSolutionDiagnosticPattynStokes(solution);
    44524451        }
    44534452}
    … …  
    44794478
    44804479        /*Get dof list: */
    4481         GetDofList(&doflist,MacAyealApproximationEnum);
     4480        GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum);
    44824481
    44834482        /*Use the dof list to index into the solution vector: */
    … …  
    45724571
    45734572        /*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
    4574         GetDofList(&doflistp,PattynApproximationEnum);
    4575         penta->GetDofList(&doflistm,MacAyealApproximationEnum);
     4573        GetDofList(&doflistp,PattynApproximationEnum,GsetEnum);
     4574        penta->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);
    45764575
    45774576        /*Get node data: */
    … …  
    46574656
    46584657        /*Get dof list: */
    4659         GetDofList(&doflist,PattynApproximationEnum);
     4658        GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    46604659
    46614660        /*Get node data: */
    … …  
    47124711}
    47134712
    4714 /*}}}*/
    4715 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{1*/
    4716 void  Penta::InputUpdateFromSolutionDiagnosticPattynStokes(double* solution){
    4717 
    4718         int i;
    4719 
    4720         const int    numdofpervertexp=2;
    4721         const int    numdofpervertexs=4;
    4722         const int    numdofp=numdofpervertexp*NUMVERTICES;
    4723         const int    numdofs=numdofpervertexs*NUMVERTICES;
    4724         int*         doflistp=NULL;
    4725         int*         doflists=NULL;
    4726         double       pattyn_values[numdofp];
    4727         double       stokes_values[numdofs];
    4728         double       vx[NUMVERTICES];
    4729         double       vy[NUMVERTICES];
    4730         double       vz[NUMVERTICES];
    4731         double       vel[NUMVERTICES];
    4732         int          dummy;
    4733         double       pressure[NUMVERTICES];
    4734         double       xyz_list[NUMVERTICES][3];
    4735 
    4736         double *vz_ptr         = NULL;
    4737         Penta  *penta          = NULL;
    4738 
    4739         /*OK, we have to add results of this element for pattyn
    4740          * and results from the penta at base for macayeal. Now recover results*/
    4741 
    4742         /*Get dof listof this element (pattyn dofs and stokes dofs): */
    4743         GetDofList(&doflistp,PattynApproximationEnum);
    4744         GetDofList(&doflists,StokesApproximationEnum);
    4745 
    4746         /*Get node data: */
    4747         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4748 
    4749         /*Use the dof list to index into the solution vector: */
    4750         for(i=0;i<numdofp;i++){
    4751                 pattyn_values[i]=solution[doflistp[i]];
    4752         }
    4753         for(i=0;i<numdofs;i++){
    4754                 stokes_values[i]=solution[doflists[i]];
    4755         }
    4756 
    4757         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    4758         for(i=0;i<NUMVERTICES;i++){
    4759                 vx[i]=stokes_values[i*numdofpervertexs+0]+pattyn_values[i*numdofpervertexp+0];
    4760                 vy[i]=stokes_values[i*numdofpervertexs+1]+pattyn_values[i*numdofpervertexp+1];
    4761                 vz[i]=stokes_values[i*numdofpervertexs+2];
    4762                 pressure[i]=stokes_values[i*numdofpervertexs+3];
    4763         }
    4764 
    4765         /*Get Vz*/
    4766         Input* vz_input=inputs->GetInput(VzEnum);
    4767         if (vz_input){
    4768                 if (vz_input->Enum()!=PentaVertexInputEnum){
    4769                         ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumToString(vz_input->Enum()));
    4770                 }
    4771                 vz_input->GetValuesPtr(&vz_ptr,&dummy);
    4772                 for(i=0;i<NUMVERTICES;i++) vz[i]=vz[i]+vz_ptr[i];
    4773         }
    4774         else{
    4775                 ISSMERROR("Cannot compute Vel there is no Vz input");
    4776         }
    4777 
    4778         /*Now Compute vel*/
    4779         for(i=0;i<NUMVERTICES;i++) {
    4780                 vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    4781         }
    4782 
    4783         /*Now, we have to move the previous Vx and Vy inputs  to old
    4784          * status, otherwise, we'll wipe them off: */
    4785         this->inputs->ChangeEnum(VxEnum,VxOldEnum);
    4786         this->inputs->ChangeEnum(VyEnum,VyOldEnum);
    4787         this->inputs->ChangeEnum(VzEnum,VzOldEnum);
    4788         this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
    4789 
    4790         /*Add vx and vy as inputs to the tria element: */
    4791         this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
    4792         this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
    4793         this->inputs->AddInput(new PentaVertexInput(VzEnum,vz));
    4794         this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
    4795         this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
    4796 
    4797         /*Free ressources:*/
    4798         xfree((void**)&doflistp);
    4799         xfree((void**)&doflists);
    4800 }
    48014713/*}}}*/
    48024714/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHutter {{{1*/
    … …  
    48214733
    48224734        /*Get dof list: */
    4823         GetDofList(&doflist);
     4735        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    48244736
    48254737        /*Get node data: */
    … …  
    49004812
    49014813        /*Get dof list: */
    4902         GetDofList(&doflist);
     4814        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    49034815
    49044816        /*Get node data: */
    … …  
    49714883
    49724884        /*Get dof list: */
    4973         GetDofList(&doflist);
     4885        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    49744886
    49754887        /*Use the dof list to index into the solution vector: */
    … …  
    50294941
    50304942        /*Get dof list: */
    5031         GetDofList(&doflist);
     4943        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    50324944
    50334945        /*Use the dof list to index into the solution vector: */
    … …  
    50674979
    50684980        /*Get dof list: */
    5069         GetDofList(&doflist);
     4981        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    50704982
    50714983        /*Use the dof list to index into the solution vector: */
    … …  
    51025014
    51035015        /*Get dof list: */
    5104         GetDofList(&doflist);
     5016        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    51055017
    51065018        /*Use the dof list to index into the solution vector: */
    … …  
    51365048
    51375049        /*Get dof list: */
    5138         GetDofList(&doflist);
     5050        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    51395051
    51405052        /*Use the dof list to index into the solution vector: */
    … …  
    51645076
    51655077        /*Get dof list: */
    5166         GetDofList(&doflist);
     5078        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    51675079
    51685080        /*Use the dof list to index into the solution vector and extrude it */

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

    r5769 r5772  
    7575                void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    7676                double RegularizeInversion(void);
    77                 void   CreateKMatrix(Mat Kgg);
    78                 void   CreatePVector(Vec pg);
     77                void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     78                void   CreatePVector(Vec pg, Vec pf);
    7979                void   DeleteResults(void);
    8080                int    GetNodeIndex(Node* node);
    … …  
    144144                void      CreatePVectorSlope( Vec pg);
    145145                void      CreatePVectorThermal( Vec pg);
    146                 void      GetDofList(int** pdoflist,int approximation_enum=0);
     146                void      GetDofList(int** pdoflist,int approximation_enum,int setenum);
    147147                void      GetDofList1(int* doflist);
    148                 int     GetElementType(void);
    149                 void    GetParameterListOnVertices(double* pvalue,int enumtype);
    150                 void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
    151                 void    GetParameterValue(double* pvalue,Node* node,int enumtype);
     148                int       GetElementType(void);
     149                void      GetParameterListOnVertices(double* pvalue,int enumtype);
     150                void      GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
     151                void      GetParameterValue(double* pvalue,Node* node,int enumtype);
    152152                void      GetPhi(double* phi, double*  epsilon, double viscosity);
    153153                void      GetSolutionFromInputsDiagnosticHoriz(Vec solutiong);

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

    r5749 r5772  
    693693/*}}}*/
    694694/*FUNCTION Tria::CreateKMatrix {{{1*/
    695 void  Tria::CreateKMatrix(Mat Kgg){
     695void  Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
    696696
    697697        /*retrive parameters: */
    … …  
    705705        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    706706        if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){
    707                 CreateKMatrixDiagnosticMacAyeal( Kgg);
     707                CreateKMatrixDiagnosticMacAyeal( Kgg,Kff,Kfs);
    708708        }
    709709        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    739739/*}}}*/
    740740/*FUNCTION Tria::CreatePVector {{{1*/
    741 void  Tria::CreatePVector(Vec pg){
     741void  Tria::CreatePVector(Vec pg, Vec pf){
    742742
    743743        int analysis_type;
    … …  
    24402440        /* Get node coordinates and dof list: */
    24412441        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2442         GetDofList(&doflist);
     2442        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    24432443
    24442444        /*retrieve some parameters: */
    … …  
    25782578        /* Get node coordinates and dof list: */
    25792579        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2580         GetDofList(&doflist);
     2580        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    25812581        this->parameters->FindParam(&dim,DimEnum);
    25822582
    … …  
    26802680        /* Get node coordinates and dof list: */
    26812681        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2682         GetDofList(&doflist);
     2682        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    26832683
    26842684        /*Retrieve all inputs we will be needed*/
    … …  
    27942794/*}}}*/
    27952795/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/
    2796 void  Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg){
     2796void  Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs){
     2797       
     2798        const int      numdof         = 2 *NUMVERTICES;
     2799        double        *Ke_gg_viscous  = NULL;
     2800        double        *Ke_gg_friction = NULL;
     2801        ElementMatrix *Ke             = NULL;
     2802       
     2803        /*compute stiffness matrices on the g-set, at the element level: */
     2804        Ke_gg_viscous =  CreateKMatrixDiagnosticMacAyealViscous();
     2805        Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction();
     2806
     2807        /*Initialize element matrix: */
     2808        Ke=this->NewElementMatrix(MacAyealApproximationEnum);
     2809
     2810        /*Add Ke_gg values to Ke element stifness matrix: */
     2811        if(Ke_gg_viscous) Ke->AddValues(Ke_gg_viscous);
     2812        if(Ke_gg_friction)Ke->AddValues(Ke_gg_friction);
     2813
     2814        /*Add Ke element stiffness matrix to global stiffness matrix: */
     2815        Ke->AddToGlobal(Kgg,Kff,Kfs);
     2816
     2817        /*Free ressources:*/
     2818        delete Ke;
     2819        xfree((void**)&Ke_gg_viscous);
     2820        xfree((void**)&Ke_gg_friction);
     2821
     2822}
     2823
     2824/*}}}*/
     2825/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/
     2826double* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){
    27972827
    27982828        /* local declarations */
    27992829        int             i,j;
     2830       
     2831        /*output: */
     2832        double*  Ke_gg=NULL;
    28002833
    28012834        /* node data: */
    28022835        const int    numdof=2*NUMVERTICES;
    28032836        double       xyz_list[NUMVERTICES][3];
    2804         int*         doflist=NULL;
    28052837
    28062838        /* gaussian points: */
    … …  
    28272859
    28282860        /* local element matrices: */
    2829         double Ke_gg[numdof][numdof]={0.0};
    28302861        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    28312862
    … …  
    28342865        /*input parameters for structural analysis (diagnostic): */
    28352866        double  thickness;
     2867       
    28362868
    28372869        /*retrieve some parameters: */
    … …  
    28392871
    28402872        /*First, if we are on water, return empty matrix: */
    2841         if(IsOnWater()) return;
     2873        if(IsOnWater()) return Ke_gg;
     2874
     2875        /*allocate output: */
     2876        Ke_gg=(double*)xcalloc(numdof*numdof,sizeof(double));
    28422877
    28432878        /* Get node coordinates and dof list: */
    28442879        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2845         GetDofList(&doflist,MacAyealApproximationEnum);
    28462880
    28472881        /*Retrieve all inputs we will be needing: */
    … …  
    28922926
    28932927                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    2894                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    2895 
    2896         }
    2897 
    2898         /*Add Ke_gg to global matrix Kgg: */
    2899         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    2900 
    2901         /*Do not forget to include friction: */
    2902         if(!IsOnShelf()){
    2903                 CreateKMatrixDiagnosticMacAyealFriction(Kgg);
    2904         }
    2905 
     2928                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) *(Ke_gg+numdof*j+i)+=Ke_gg_gaussian[i][j];
     2929
     2930        }
    29062931
    29072932        /*Clean up and return*/
    29082933        delete gauss;
    2909         xfree((void**)&doflist);
    2910 }
    2911 /*}}}*/
    2912 /*FUNCTION Tria::CreateKMatrixCouplingMacAyealPattynFriction {{{1*/
    2913 void  Tria::CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg){
     2934        return Ke_gg;
     2935}
     2936/*}}}*/
     2937/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
     2938void  Tria::CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs){
     2939       
     2940        const int      numdof         = 2 *NUMVERTICES;
     2941        double        *Ke_gg_friction = NULL;
     2942        ElementMatrix *Ke             = NULL;
     2943       
     2944        /*compute stiffness matrices on the g-set, at the element level: */
     2945        Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction();
     2946
     2947        if(Ke_gg_friction){
     2948
     2949                /*Initialize element matrix: */
     2950                Ke=this->NewElementMatrix(MacAyealApproximationEnum);
     2951
     2952                /*Add Ke_gg values to Ke element stifness matrix: */
     2953                Ke->AddValues(Ke_gg_friction);
     2954       
     2955                /*Add Ke element stiffness matrix to global stiffness matrix: */
     2956                Ke->AddToGlobal(Kgg,Kff,Kfs);
     2957       
     2958                /*Free ressources:*/
     2959                delete Ke;
     2960        }
     2961
     2962        /*Free ressources:*/
     2963        xfree((void**)&Ke_gg_friction);
     2964}
     2965
     2966/*}}}*/
     2967/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
     2968double*  Tria::CreateKMatrixDiagnosticMacAyealFriction(void){
    29142969
    29152970        /* local declarations */
    … …  
    29172972        int       analysis_type;
    29182973
     2974        /*output: */
     2975        double*   Ke_gg=NULL;
     2976
    29192977        /* node data: */
    29202978        const int numdof   = 2 *NUMVERTICES;
    29212979        double    xyz_list[NUMVERTICES][3];
    2922         int*      doflistm=NULL;
    2923         int*      doflistp=NULL;
    29242980        int       numberofdofspernode=2;
    29252981
    … …  
    29342990
    29352991        /* local element matrices: */
    2936         double Ke_gg[numdof][numdof]={0.0};
    29372992        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    29382993       
    … …  
    29523007        /*inputs: */
    29533008        int  drag_type;
     3009       
    29543010
    29553011        /*retrive parameters: */
    … …  
    29653021        /* Get node coordinates and dof list: */
    29663022        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2967         GetDofList(&doflistm,MacAyealApproximationEnum);
    2968         GetDofList(&doflistp,PattynApproximationEnum);
    29693023
    29703024        if (IsOnShelf()){
    29713025                /*no friction, do nothing*/
    2972                 return;
    2973         }
     3026                return Ke_gg;
     3027        }
     3028
     3029        /*allocte Ke_gg matrix: */
     3030        Ke_gg=(double*)xcalloc(numdof*numdof,sizeof(double));
    29743031
    29753032        /*build friction object, used later on: */
    … …  
    30133070                                        &Ke_gg_gaussian[0][0],0);
    30143071
    3015                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    3016 
    3017         }
    3018 
    3019         /*Add Ke_gg to global matrix Kgg: */
    3020         MatSetValues(Kgg,numdof,doflistm,numdof,doflistp,(const double*)Ke_gg,ADD_VALUES);
    3021         MatSetValues(Kgg,numdof,doflistp,numdof,doflistm,(const double*)Ke_gg,ADD_VALUES);
     3072                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) *(Ke_gg+numdof*i+j)+=Ke_gg_gaussian[i][j];
     3073
     3074        }
    30223075
    30233076        /*Clean up and return*/
    30243077        delete gauss;
    30253078        delete friction;
    3026         xfree((void**)&doflistm);
    3027         xfree((void**)&doflistp);
    3028 }       
    3029 /*}}}*/
    3030 /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
    3031 void  Tria::CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg){
     3079        return Ke_gg;
     3080}
     3081/*}}}*/
     3082/*FUNCTION Tria::CreateKMatrixCouplingMacAyealPattynFriction {{{1*/
     3083void  Tria::CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg){
    30323084
    30333085        /* local declarations */
    … …  
    30383090        const int numdof   = 2 *NUMVERTICES;
    30393091        double    xyz_list[NUMVERTICES][3];
    3040         int*      doflist=NULL;
     3092        int*      doflistm=NULL;
     3093        int*      doflistp=NULL;
    30413094        int       numberofdofspernode=2;
    30423095
    … …  
    30823135        /* Get node coordinates and dof list: */
    30833136        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3084         GetDofList(&doflist,MacAyealApproximationEnum);
     3137        GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);
     3138        GetDofList(&doflistp,PattynApproximationEnum,GsetEnum);
    30853139
    30863140        if (IsOnShelf()){
    … …  
    31343188
    31353189        /*Add Ke_gg to global matrix Kgg: */
    3136         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     3190        MatSetValues(Kgg,numdof,doflistm,numdof,doflistp,(const double*)Ke_gg,ADD_VALUES);
     3191        MatSetValues(Kgg,numdof,doflistp,numdof,doflistm,(const double*)Ke_gg,ADD_VALUES);
    31373192
    31383193        /*Clean up and return*/
    31393194        delete gauss;
    31403195        delete friction;
    3141         xfree((void**)&doflist);
     3196        xfree((void**)&doflistm);
     3197        xfree((void**)&doflistp);
    31423198}       
    31433199/*}}}*/
    … …  
    31963252        /* Get node coordinates and dof list: */
    31973253        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3198         GetDofList(&doflist,PattynApproximationEnum);
     3254        GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    31993255
    32003256        if (IsOnShelf()){
    … …  
    32693325        if(IsOnWater())return;
    32703326
    3271         GetDofList(&doflist);
     3327        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    32723328
    32733329        /*Spawn 3 sing elements: */
    … …  
    33193375        /* Get node coordinates and dof list: */
    33203376        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3321         GetDofList(&doflist);
     3377        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    33223378
    33233379        /*Build normal vector to the surface:*/
    … …  
    34183474        /* Get node coordinates and dof list: */
    34193475        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3420         GetDofList(&doflist);
     3476        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    34213477
    34223478
    … …  
    35073563        /* Get node coordinates and dof list: */
    35083564        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3509         GetDofList(&doflist);
     3565        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    35103566
    35113567        /*Retrieve all inputs we will be needing: */
    … …  
    36663722        /* Get node coordinates and dof list: */
    36673723        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3668         GetDofList(&doflist);
     3724        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    36693725
    36703726        /*Retrieve all inputs we will be needed*/
    … …  
    37503806
    37513807        GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes);
    3752         GetDofList(&doflist);
     3808        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    37533809
    37543810        /* Start looping on the number of gaussian points: */
    … …  
    38143870        /* Get node coordinates and dof list: */
    38153871        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3816         GetDofList(&doflist);
     3872        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    38173873
    38183874        //recover material parameters
    … …  
    38893945        /* Get node coordinates and dof list: */
    38903946        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3891         GetDofList(&doflist);
     3947        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    38923948
    38933949        /*retrieve inputs :*/
    … …  
    39544010        /* Get node coordinates and dof list: */
    39554011        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3956         GetDofList(&doflist);
     4012        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    39574013
    39584014        /*Retrieve all inputs we will be needing: */
    … …  
    40184074        /* Get node coordinates and dof list: */
    40194075        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4020         GetDofList(&doflist);
     4076        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    40214077
    40224078        /*Retrieve all inputs we will be needing: */
    … …  
    40884144        /* Get node coordinates and dof list: */
    40894145        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4090         GetDofList(&doflist);
     4146        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    40914147
    40924148        /*Retrieve all inputs we will be needing: */
    … …  
    41734229        /* Get node coordinates and dof list: */
    41744230        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4175         GetDofList(&doflist,MacAyealApproximationEnum);
     4231        GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum);
    41764232
    41774233        /* Start  looping on the number of gaussian points: */
    … …  
    42574313        /* Get node coordinates and dof list: */
    42584314        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4259         GetDofList(&doflist);
     4315        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    42604316
    42614317        /*Retrieve all inputs we will be needing: */
    … …  
    43504406        /* Get node coordinates and dof list: */
    43514407        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4352         GetDofList(&doflist);
     4408        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    43534409
    43544410        /* Recover input data: */
    … …  
    45714627        /* Get node coordinates and dof list: */
    45724628        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4573         GetDofList(&doflist);
     4629        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    45744630
    45754631        /* Recover input data: */
    … …  
    47554811        if(IsOnWater())return;
    47564812
    4757         GetDofList(&doflist);
     4813        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    47584814       
    47594815        /* Get parameters */
    … …  
    48324888        /* Get node coordinates and dof list: */
    48334889        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4834         GetDofList(&doflist);
     4890        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    48354891
    48364892        /*Retrieve all inputs we will be needing: */
    … …  
    49024958        /* Get node coordinates and dof list: */
    49034959        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4904         GetDofList(&doflist);
     4960        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    49054961
    49064962        /*Retrieve all inputs we will be needing: */
    … …  
    49735029        /* Get node coordinates and dof list: */
    49745030        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4975         GetDofList(&doflist);
     5031        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    49765032
    49775033        /*Retrieve all inputs we will be needing: */
    … …  
    50525108        /* Get node coordinates and dof list: */
    50535109        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    5054         GetDofList(&doflist);
     5110        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    50555111
    50565112        //recover material parameters
    … …  
    51455201        /* Get node coordinates and dof list: */
    51465202        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    5147         GetDofList(&doflist);
     5203        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    51485204
    51495205        //recover material parameters
    … …  
    52315287/*}}}*/
    52325288/*FUNCTION Tria::GetDofList {{{1*/
    5233 void  Tria::GetDofList(int** pdoflist, int approximation_enum){
     5289void  Tria::GetDofList(int** pdoflist, int approximation_enum,int setenum){
    52345290
    52355291        int i,j;
    … …  
    52425298        /*First, figure out size of doflist: */
    52435299        for(i=0;i<3;i++){
    5244                 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum);
     5300                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    52455301        }
    52465302
    … …  
    52515307        count=0;
    52525308        for(i=0;i<3;i++){
    5253                 nodes[i]->GetDofList(doflist+count,approximation_enum);
    5254                 count+=nodes[i]->GetNumberOfDofs(approximation_enum);
     5309                nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
     5310                count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    52555311        }
    52565312
    … …  
    52585314        *pdoflist=doflist;
    52595315
     5316}
     5317/*}}}*/
     5318/*FUNCTION Tria::GetInternalDofList {{{1*/
     5319int*  Tria::GetInternalDofList(int approximation_enum,int setenum){
     5320
     5321        int i,j;
     5322        int numberofdofs=0;
     5323        int count=0;
     5324        int count2=0;
     5325        int ndof;
     5326        int ngdof;
     5327
     5328        /*output: */
     5329        int* doflist=NULL;
     5330
     5331        /*First, figure out size of doflist: */
     5332        for(i=0;i<3;i++){
     5333                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     5334        }
     5335
     5336        if(numberofdofs){
     5337                /*Allocate: */
     5338                doflist=(int*)xmalloc(numberofdofs*sizeof(int));
     5339
     5340                /*Populate: */
     5341                count=0;
     5342                for(i=0;i<3;i++){
     5343                        nodes[i]->GetInternalDofList(doflist+count,approximation_enum,setenum);
     5344                        count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     5345                }
     5346               
     5347                /*We now have something like: [0 1 0 2 0 1]. Offset by gsize, to get something like:
     5348                 * [0 1 2 4 5 6]:*/
     5349                count=0;
     5350                count2=0;
     5351                for(i=0;i<3;i++){
     5352                        ndof=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     5353                        for(j=count;j<count+ndof;j++)doflist[j]+=count2;
     5354                        count+=ndof;
     5355                        count2+=nodes[i]->GetNumberOfDofs(approximation_enum,GsetEnum);
     5356                }
     5357
     5358        }
     5359        else doflist=NULL;
     5360
     5361        return doflist;
     5362
     5363}
     5364/*}}}*/
     5365/*FUNCTION Tria::GetExternalDofList {{{1*/
     5366int* Tria::GetExternalDofList(int approximation_enum,int setenum){
     5367
     5368        int i,j;
     5369        int numberofdofs=0;
     5370        int count=0;
     5371
     5372        /*output: */
     5373        int* doflist=NULL;
     5374
     5375        /*First, figure out size of doflist: */
     5376        for(i=0;i<3;i++){
     5377                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     5378        }
     5379
     5380        if(numberofdofs){
     5381                /*Allocate: */
     5382                doflist=(int*)xmalloc(numberofdofs*sizeof(int));
     5383
     5384                /*Populate: */
     5385                count=0;
     5386                for(i=0;i<3;i++){
     5387                        nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
     5388                        count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     5389                }
     5390        }
     5391        else doflist=NULL;
     5392
     5393        return doflist;
     5394}
     5395/*}}}*/
     5396/*FUNCTION Tria::GetNumberOfDofs {{{1*/
     5397int Tria::GetNumberOfDofs(int approximation_enum,int setenum){
     5398
     5399        int i;
     5400       
     5401        /*output: */
     5402        int numberofdofs=0;
     5403
     5404        for(i=0;i<3;i++){
     5405                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     5406        }
     5407
     5408        return numberofdofs;
    52605409}
    52615410/*}}}*/
    … …  
    53505499
    53515500        /*Get dof list: */
    5352         GetDofList(&doflist);
     5501        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    53535502
    53545503        /*Get inputs*/
    … …  
    53935542
    53945543        /*Get dof list: */
    5395         GetDofList(&doflist);
     5544        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    53965545
    53975546        /*Get inputs*/
    … …  
    55865735
    55875736        /*Get dof list: */
    5588         GetDofList(&doflist);
     5737        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    55895738
    55905739        /*Use the dof list to index into the solution vector: */
    … …  
    56195768
    56205769        /*Get dof list: */
    5621         GetDofList(&doflist);
     5770        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    56225771
    56235772        /*Use the dof list to index into the solution vector: */
    … …  
    56605809       
    56615810        /*Get dof list: */
    5662         GetDofList(&doflist);
     5811        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    56635812
    56645813        /*Use the dof list to index into the solution vector: */
    … …  
    57245873       
    57255874        /*Get dof list: */
    5726         GetDofList(&doflist);
     5875        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    57275876
    57285877        /*Use the dof list to index into the solution vector: */
    … …  
    57895938
    57905939        /*Get dof list: */
    5791         GetDofList(&doflist);
     5940        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    57925941
    57935942        /*Use the dof list to index into the solution vector: */
    … …  
    58305979}
    58315980/*}}}*/
     5981/*FUNCTION Tria::NewElementMatrix{{{1*/
     5982ElementMatrix* Tria::NewElementMatrix(int approximation){
     5983
     5984        /*parameters: */
     5985        bool kff=false;
     5986
     5987        /*output: */
     5988        ElementMatrix* Ke=NULL;
     5989        int gsize;
     5990        int fsize;
     5991        int ssize;
     5992        int* gexternaldoflist=NULL;
     5993        int* finternaldoflist=NULL;
     5994        int* fexternaldoflist=NULL;
     5995        int* sinternaldoflist=NULL;
     5996        int* sexternaldoflist=NULL;
     5997        bool symmetric=true;
     5998
     5999        /*retrieve some parameters: */
     6000        this->parameters->FindParam(&kff,KffEnum);
     6001
     6002        /*get number of dofs in sets g,f and s: */
     6003        gsize=this->GetNumberOfDofs(approximation,GsetEnum);
     6004        if(kff){
     6005                fsize=this->GetNumberOfDofs(approximation,FsetEnum);
     6006                ssize=this->GetNumberOfDofs(approximation,SsetEnum);
     6007        }
     6008
     6009        /*get dof lists for f and s set: */
     6010        if(!kff) gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum);
     6011        else{
     6012                finternaldoflist=this->GetInternalDofList(approximation,FsetEnum);
     6013                fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum);
     6014                sinternaldoflist=this->GetInternalDofList(approximation,SsetEnum);
     6015                sexternaldoflist=this->GetExternalDofList(approximation,SsetEnum);
     6016        }
     6017
     6018        /*Use symmetric constructor for ElementMatrix: */
     6019        if(!kff) Ke=new ElementMatrix(gsize,symmetric,gexternaldoflist);
     6020        else     Ke=new ElementMatrix(gsize,symmetric,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize);
     6021
     6022        /*Free ressources and return:*/
     6023        xfree((void**)&gexternaldoflist);
     6024        xfree((void**)&finternaldoflist);
     6025        xfree((void**)&fexternaldoflist);
     6026        xfree((void**)&sinternaldoflist);
     6027        xfree((void**)&sexternaldoflist);
     6028
     6029        return Ke;
     6030}
     6031/*}}}*/
     6032/*FUNCTION Tria::NewElementVector{{{1*/
     6033ElementVector* Tria::NewElementVector(int approximation){
     6034
     6035        /*parameters: */
     6036        bool kff=false;
     6037
     6038        /*output: */
     6039        ElementVector* pe=NULL;
     6040        int gsize;
     6041        int fsize;
     6042        int* gexternaldoflist=NULL;
     6043        int* finternaldoflist=NULL;
     6044        int* fexternaldoflist=NULL;
     6045
     6046        /*retrieve some parameters: */
     6047        this->parameters->FindParam(&kff,KffEnum);
     6048
     6049        /*get number of dofs in sets g,f and s: */
     6050        gsize=this->GetNumberOfDofs(approximation,GsetEnum);
     6051        if(kff)fsize=this->GetNumberOfDofs(approximation,FsetEnum);
     6052
     6053        /*get dof lists for f and s set: */
     6054        if(!kff){
     6055                gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum);
     6056        }
     6057        else{
     6058                finternaldoflist=this->GetInternalDofList(approximation,FsetEnum);
     6059                fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum);
     6060        }
     6061
     6062        /*Use symmetric constructor for ElementVector: */
     6063        if(!kff)pe=new ElementVector(gsize,gexternaldoflist);
     6064        else    pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize);
     6065
     6066        /*Free ressources and return:*/
     6067        xfree((void**)&gexternaldoflist);
     6068        xfree((void**)&finternaldoflist);
     6069        xfree((void**)&fexternaldoflist);
     6070       
     6071        return pe;
     6072}
     6073/*}}}*/
    58326074/*FUNCTION Tria::SetClone {{{1*/
    58336075void  Tria::SetClone(int* minranks){

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

    r5746 r5772  
    1717class Matice;
    1818class Matpar;
     19class ElementMatrix;
     20class ElementVector;
    1921
    2022#include "../../include/include.h"
    … …  
    7274                void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    7375                double RegularizeInversion(void);
    74                 void   CreateKMatrix(Mat Kgg);
    75                 void   CreatePVector(Vec pg);
     76                void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     77                void   CreatePVector(Vec pg, Vec pf);
    7678                int    GetNodeIndex(Node* node);
    7779                bool   IsOnBed();
    … …  
    122124                void      CreateKMatrixBalancedthickness_CG(Mat Kgg);
    123125                void      CreateKMatrixBalancedvelocities(Mat Kgg);
    124                 void      CreateKMatrixDiagnosticMacAyeal(Mat Kgg);
     126                void      CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs);
     127                double*   CreateKMatrixDiagnosticMacAyealViscous(void);
     128                double*   CreateKMatrixDiagnosticMacAyealFriction(void);
     129                void      CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs);
    125130                void      CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg);
    126                 void      CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg);
    127131                void      CreateKMatrixDiagnosticPattynFriction(Mat Kgg);
    128132                void      CreateKMatrixDiagnosticHutter(Mat Kgg);
    … …  
    147151                void      CreatePVectorThermalSheet( Vec pg);
    148152                void      CreatePVectorThermalShelf( Vec pg);
    149                 double  GetArea(void);
    150                 int     GetElementType(void);
    151                 void      GetDofList(int** pdoflist,int approximation_enum=0);
     153                double    GetArea(void);
     154                int       GetElementType(void);
     155                void      GetDofList(int** pdoflist,int approximation_enum,int setenum);
    152156                void      GetDofList1(int* doflist);
    153                 void    GetParameterListOnVertices(double* pvalue,int enumtype);
    154                 void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
    155                 void    GetParameterValue(double* pvalue,Node* node,int enumtype);
     157                int*      GetInternalDofList(int approximation_enum,int setenum);
     158                int*      GetExternalDofList(int approximation_enum,int setenum);
     159                int       GetNumberOfDofs(int approximation_enum,int setenum);
     160                void      GetParameterListOnVertices(double* pvalue,int enumtype);
     161                void      GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
     162                void      GetParameterValue(double* pvalue,Node* node,int enumtype);
    156163                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    157164                void      GetSolutionFromInputsDiagnosticHutter(Vec solution);
    158                 void    GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
     165                void      GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
    159166                void      GradjDragStokes(Vec gradient);
    160167                void      InputUpdateFromSolutionAdjointBalancedthickness( double* solution);
    … …  
    164171                void      InputUpdateFromSolutionOneDof(double* solution,int enum_type);
    165172                bool      IsInput(int name);
     173                ElementMatrix* NewElementMatrix(int approximation);
     174                ElementVector* NewElementVector(int approximation);
    166175                void      SetClone(int* minranks);
    167176                void      SurfaceNormal(double* surface_normal, double xyz_list[3][3]);

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

    r5103 r5772  
    3838        analysis_type_list=(int*)xmalloc(nummodels*sizeof(int));
    3939        m_nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*));
    40         m_yg=(Vec*)xmalloc(nummodels*sizeof(Vec));
    4140        m_ys=(Vec*)xmalloc(nummodels*sizeof(Vec));
    4241
    … …  
    4443        for(i=0;i<nummodels;i++)analysis_type_list[i]=analyses[i];
    4544        for(i=0;i<nummodels;i++)m_nodesets[i]=NULL;
    46         for(i=0;i<nummodels;i++)m_yg[i]=NULL;
    4745        for(i=0;i<nummodels;i++)m_ys[i]=NULL;
    4846
    … …  
    5755                this->SetCurrentConfiguration(analysis_type);
    5856       
    59                 _printf_("      create degrees of freedom\n");
     57                _printf_("      create vertex degrees of freedom\n");
    6058                VerticesDofx(&partition,&tpartition,vertices,parameters);
     59
     60                _printf_("      resolve node constraints\n");
     61                SpcNodesx(nodes,constraints,analysis_type);
     62       
     63                _printf_("      create nodal degrees of freedom\n");
    6164                NodesDofx(nodes,parameters,analysis_type);
    62 
    63                 _printf_("      create single point constraints\n");
    64                 SpcNodesx( &m_yg[i], nodes,constraints,analysis_type);
     65       
     66                _printf_("      create nodal constraints vector\n");
     67                CreateNodalConstraintsx(&m_ys[i],nodes,analysis_type);
    6568
    6669                _printf_("      create node sets\n");
    6770                BuildNodeSetsx(&m_nodesets[i], nodes,analysis_type);
    68 
    69                 _printf_("      reducing single point constraints vector\n");
    70                 Reducevectorgtosx(&m_ys[i], m_yg[i],m_nodesets[i]);
    7171
    7272                _printf_("      configuring element and loads\n");
    … …  
    9797                NodeSets* temp_nodesets=m_nodesets[i];
    9898                delete temp_nodesets;
    99                 Vec temp_yg=m_yg[i];
    100                 VecFree(&temp_yg);
    10199                Vec temp_ys=m_ys[i];
    102100                VecFree(&temp_ys);
    … …  
    105103        /*Delete dynamically allocated arrays: */
    106104        xfree((void**)&m_nodesets);
    107         xfree((void**)&m_yg);
    108105        xfree((void**)&m_ys);
    109106
    … …  
    147144        /*activate matrices/vectors: */
    148145        nodesets=m_nodesets[analysis_counter];
    149         yg=m_yg[analysis_counter];
    150146        ys=m_ys[analysis_counter];
    151147

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

    r5057 r5772  
    4141                //multiple  sets of matrices/vectors for each analysis_type. m stands for multiple
    4242                NodeSets**           m_nodesets; //boundary conditions dof sets
    43                 Vec*                 m_yg; //boundary conditions in global g-set
    4443                Vec*                 m_ys; //boundary conditions, in reduced s-set
    4544

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

    r5631 r5772  
    188188        IoModelFetchData(&this->stokesreconditioning,iomodel_handle,"stokesreconditioning");
    189189        IoModelFetchData(&this->waitonlock,iomodel_handle,"waitonlock");
     190        IoModelFetchData(&this->kff,iomodel_handle,"kff");
    190191
    191192        /*!Get thermal parameters: */
    … …  
    324325        this->connectivity=0;
    325326        this->lowmem=0;
     327        this->kff=0;
    326328        this->optscal=NULL;
    327329        this->yts=0;

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

    r5631 r5772  
    150150                double  yts;
    151151                double  waitonlock;
     152                int     kff;
    152153
    153154                /*thermal parameters: */

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5743 r5772  
    307307/*}}}*/
    308308/*FUNCTION Icefront::CreateKMatrix {{{1*/
    309 void  Icefront::CreateKMatrix(Mat Kgg){
     309void  Icefront::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
    310310
    311311        /*No stiffness loads applied, do nothing: */
    … …  
    315315/*}}}*/
    316316/*FUNCTION Icefront::CreatePVector {{{1*/
    317 void  Icefront::CreatePVector(Vec pg){
     317void  Icefront::CreatePVector(Vec pg, Vec pf){
    318318
    319319        /*Checks in debugging mode*/
    … …  
    331331                inputs->GetParameterValue(&type,TypeEnum);
    332332                if (type==MacAyeal2dIceFrontEnum){
    333                         CreatePVectorDiagnosticMacAyeal2d( pg);
     333                        CreatePVectorDiagnosticMacAyeal2d( pg,pf);
    334334                }
    335335                else if (type==MacAyeal3dIceFrontEnum){
    … …  
    353353/*}}}*/
    354354/*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/
    355 void  Icefront::PenaltyCreateKMatrix(Mat Kgg,double kmax){
     355void  Icefront::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs, double kmax){
    356356        /*do nothing: */
    357357}
    358358/*}}}*/
    359359/*FUNCTION Icefront::PenaltyCreatePVector{{{1*/
    360 void  Icefront::PenaltyCreatePVector(Vec pg,double kmax){
     360void  Icefront::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){
    361361        /*do nothing: */
    362362}
    … …  
    423423/*Icefront numerics: */
    424424/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{1*/
    425 void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg){
     425void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg,Vec pf){
    426426
    427427        /*Constants*/
    428428        const int numnodes= NUMVERTICESSEG;
    429429        const int numdofs = numnodes *NDOF2;
     430
     431        /*output: */
     432        ElementVector* pe=NULL;
    430433
    431434        /*Intermediary*/
    … …  
    434437        double     thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity;
    435438        double     water_pressure,air_pressure,surface_under_water,base_under_water;
    436         int       *doflist = NULL;
    437439        double     xyz_list[numnodes][3];
    438440        double     pe_g[numdofs] = {0.0};
    … …  
    444446
    445447        /* Get dof list and node coordinates: */
    446         GetDofList(&doflist,MacAyealApproximationEnum);
    447448        GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
    448449
    … …  
    493494                }
    494495        }
    495        
    496         VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
     496
     497        /*Initialize element matrix: */
     498        pe=this->NewElementVector(MacAyealApproximationEnum);
     499
     500        /*Add pe_g values to pe element stifness load: */
     501        pe->AddValues(&pe_g[0]);
     502
     503        /*Add pe element load vector to global load vector: */
     504        pe->AddToGlobal(pg,pf);
     505
     506        /*Free ressources:*/
     507        delete pe;
    497508
    498509        /*Free ressources:*/
    499510        delete gauss;
    500         xfree((void**)&doflist);
    501511}
    502512/*}}}*/
    … …  
    521531        icefront->element=tria;
    522532        icefront->inputs->AddInput(new IntInput(TypeEnum,MacAyeal2dIceFrontEnum));
    523         icefront->CreatePVectorDiagnosticMacAyeal2d(pg);
     533        icefront->CreatePVectorDiagnosticMacAyeal2d(pg,NULL);
    524534
    525535        /*clean-up*/
    … …  
    574584
    575585        /* Get dof list and node coordinates: */
    576         GetDofList(&doflist,PattynApproximationEnum);
     586        GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    577587        GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes);
    578588       
    … …  
    692702
    693703        /* Get dof list and node coordinates: */
    694         GetDofList(&doflist,StokesApproximationEnum);
     704        GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
    695705        GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes);
    696706       
    … …  
    766776/*}}}*/
    767777/*FUNCTION Icefront::GetDofList {{{1*/
    768 void  Icefront::GetDofList(int** pdoflist,int approximation_enum){
     778void  Icefront::GetDofList(int** pdoflist,int approximation_enum,int setenum){
    769779
    770780        int i,j;
    … …  
    792802        /*Figure out size of doflist: */
    793803        for(i=0;i<numberofnodes;i++){
    794                 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum);
     804                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    795805        }
    796806
    … …  
    801811        count=0;
    802812        for(i=0;i<numberofnodes;i++){
    803                 nodes[i]->GetDofList(doflist+count,approximation_enum);
    804                 count+=nodes[i]->GetNumberOfDofs(approximation_enum);
     813                nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
     814                count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    805815        }
    806816
    807817        /*Assign output pointers:*/
    808818        *pdoflist=doflist;
     819}
     820/*}}}*/
     821/*FUNCTION Icefront::GetInternalDofList {{{1*/
     822int* Icefront::GetInternalDofList(int approximation_enum,int setenum){
     823
     824        int i,j;
     825        int numberofdofs=0;
     826        int count=0;
     827        int count2=0;
     828        int type;
     829        int numberofnodes=2;
     830        int ndof;
     831
     832        /*output: */
     833        int* doflist=NULL;
     834
     835        /*recover type: */
     836        inputs->GetParameterValue(&type,TypeEnum);
     837
     838        /*Some checks for debugging*/
     839        ISSMASSERT(nodes);
     840               
     841        /*How many nodes? :*/
     842        if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
     843         numberofnodes=2;
     844        else
     845         numberofnodes=4;
     846       
     847        /*Figure out size of doflist: */
     848        for(i=0;i<numberofnodes;i++){
     849                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     850        }
     851
     852        if(numberofdofs){
     853                /*Allocate: */
     854                doflist=(int*)xmalloc(numberofdofs*sizeof(int));
     855
     856                /*Populate: */
     857                count=0;
     858                for(i=0;i<numberofnodes;i++){
     859                        nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
     860                        count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     861                }
     862               
     863                /*We now have something like: [0 1 0 1]. Offset by gsize, to get something like:
     864                 * [0 1 2 3]: */
     865                count=0;
     866                count2=0;
     867                for(i=0;i<numberofnodes;i++){
     868                        ndof=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     869                        for(j=count;j<count+ndof;j++)doflist[j]+=count2;
     870                        count+=ndof;
     871                        count2+=nodes[i]->GetNumberOfDofs(approximation_enum,GsetEnum);
     872                }
     873        }
     874
     875        return doflist;
     876}
     877/*}}}*/
     878/*FUNCTION Icefront::GetExternalDofList{{{1*/
     879int*  Icefront::GetExternalDofList(int approximation_enum,int setenum){
     880
     881        int i,j;
     882        int numberofdofs=0;
     883        int count=0;
     884        int type;
     885        int numberofnodes=2;
     886
     887        /*output: */
     888        int* doflist=NULL;
     889
     890        /*recover type: */
     891        inputs->GetParameterValue(&type,TypeEnum);
     892
     893        /*Some checks for debugging*/
     894        ISSMASSERT(nodes);
     895               
     896        /*How many nodes? :*/
     897        if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
     898         numberofnodes=2;
     899        else
     900         numberofnodes=4;
     901       
     902        /*Figure out size of doflist: */
     903        for(i=0;i<numberofnodes;i++){
     904                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     905        }
     906
     907        if(numberofdofs){
     908                /*Allocate: */
     909                doflist=(int*)xmalloc(numberofdofs*sizeof(int));
     910
     911                /*Populate: */
     912                count=0;
     913                for(i=0;i<numberofnodes;i++){
     914                        nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
     915                        count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     916                }
     917        }
     918
     919        return doflist;
     920}
     921/*}}}*/
     922/*FUNCTION Icefront::GetNumberOfDofs {{{1*/
     923int Icefront::GetNumberOfDofs(int approximation_enum,int setenum){
     924       
     925       
     926        int i;
     927        int numberofdofs=0;
     928        int type;
     929        int numberofnodes=2;
     930
     931        /*recover type: */
     932        inputs->GetParameterValue(&type,TypeEnum);
     933
     934        /*Some checks for debugging*/
     935        ISSMASSERT(nodes);
     936               
     937        /*How many nodes? :*/
     938        if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
     939         numberofnodes=2;
     940        else
     941         numberofnodes=4;
     942       
     943        /*Figure out size of doflist: */
     944        for(i=0;i<numberofnodes;i++){
     945                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     946        }
     947
     948        return numberofdofs;
    809949}
    810950/*}}}*/
    … …  
    13091449}
    13101450/*}}}*/
     1451/*FUNCTION Icefront::NewElementVector{{{1*/
     1452ElementVector* Icefront::NewElementVector(int approximation){
     1453
     1454        /*parameters: */
     1455        bool kff=false;
     1456
     1457        /*output: */
     1458        ElementVector* pe=NULL;
     1459        int gsize;
     1460        int fsize;
     1461        int* gexternaldoflist=NULL;
     1462        int* finternaldoflist=NULL;
     1463        int* fexternaldoflist=NULL;
     1464
     1465        /*retrieve some parameters: */
     1466        this->parameters->FindParam(&kff,KffEnum);
     1467
     1468        /*get number of dofs in sets g and f: */
     1469        gsize=this->GetNumberOfDofs(approximation,GsetEnum);
     1470        if(kff)fsize=this->GetNumberOfDofs(approximation,FsetEnum);
     1471
     1472        /*get dof lists for f and s set: */
     1473        if(!kff){
     1474                gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum);
     1475        }
     1476        else{
     1477                finternaldoflist=this->GetInternalDofList(approximation,FsetEnum);
     1478                fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum);
     1479        }
     1480
     1481        /*Use symmetric constructor for ElementVector: */
     1482        if(!kff)pe=new ElementVector(gsize,gexternaldoflist);
     1483        else    pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize);
     1484
     1485        /*Free ressources and return:*/
     1486        xfree((void**)&gexternaldoflist);
     1487        xfree((void**)&finternaldoflist);
     1488        xfree((void**)&fexternaldoflist);
     1489       
     1490        return pe;
     1491}
     1492/*}}}*/

  • issm/trunk/src/c/objects/Loads/Icefront.h

    r5722 r5772  
    1616class Element;
    1717class IoModel;
     18class ElementVector;
    1819/*}}}*/
    1920
    … …  
    6970                void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    7071                void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    71                 void  CreateKMatrix(Mat Kgg);
    72                 void  CreatePVector(Vec pg);
    73                 void  PenaltyCreateKMatrix(Mat Kgg,double kmax);
    74                 void  PenaltyCreatePVector(Vec pg,double kmax);
     72                void  CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     73                void  CreatePVector(Vec pg, Vec pf);
     74                void  PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax);
     75                void  PenaltyCreatePVector(Vec pg,Vec pf, double kmax);
    7576                bool  InAnalysis(int analysis_type);
    7677                /*}}}*/
    7778                /*Load management: {{{1*/
    78                 void  CreatePVectorDiagnosticMacAyeal2d(Vec pg);
     79                void  CreatePVectorDiagnosticMacAyeal2d(Vec pg,Vec pf);
    7980                void  CreatePVectorDiagnosticMacAyeal3d(Vec pg);
    8081                void  CreatePVectorDiagnosticPattyn(Vec pg);
    8182                void  CreatePVectorDiagnosticStokes(Vec pg);
    82                 void  GetDofList(int** pdoflist,int approximation_enum=0);
     83                void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
     84                int*  GetInternalDofList(int approximation_enum,int setenum);
     85                int*  GetExternalDofList(int approximation_enum,int setenum);
     86                int   GetNumberOfDofs(int approximation_enum,int setenum);
    8387                void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    8488                                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
    … …  
    8690                                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
    8791                void GetNormal(double* normal,double xyz_list[2][3]);
     92                ElementVector* NewElementVector(int approximation);
    8893                /*}}}*/
    8994};

  • issm/trunk/src/c/objects/Loads/Load.h

    r4575 r5772  
    2626                virtual void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
    2727                virtual void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
    28                 virtual void  CreateKMatrix(Mat Kgg)=0;
    29                 virtual void  CreatePVector(Vec pg)=0;
    30                 virtual void  PenaltyCreateKMatrix(Mat Kgg,double kmax)=0;
    31                 virtual void  PenaltyCreatePVector(Vec pg,double kmax)=0;
     28                virtual void  CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs)=0;
     29                virtual void  CreatePVector(Vec pg, Vec pf)=0;
     30                virtual void  PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs, double kmax)=0;
     31                virtual void  PenaltyCreatePVector(Vec pg, Vec pf, double kmax)=0;
    3232                virtual bool  InAnalysis(int analysis_type)=0;
    3333                /*}}}*/

  • issm/trunk/src/c/objects/Loads/Numericalflux.cpp

    r5741 r5772  
    328328/*}}}*/
    329329/*FUNCTION Numericalflux::CreateKMatrix {{{1*/
    330 void  Numericalflux::CreateKMatrix(Mat Kgg){
     330void  Numericalflux::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){
    331331
    332332        int type;
    … …  
    349349/*}}}*/
    350350/*FUNCTION Numericalflux::CreatePVector {{{1*/
    351 void  Numericalflux::CreatePVector(Vec pg){
     351void  Numericalflux::CreatePVector(Vec pg,Vec pf){
    352352
    353353        int type;
    … …  
    374374/*}}}*/
    375375/*FUNCTION Numericalflux::PenaltyCreateKMatrix {{{1*/
    376 void  Numericalflux::PenaltyCreateKMatrix(Mat Kgg,double kmax){
     376void  Numericalflux::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){
    377377
    378378        /*No stiffness loads applied, do nothing: */
    … …  
    382382/*}}}*/
    383383/*FUNCTION Numericalflux::PenaltyCreatePVector{{{1*/
    384 void  Numericalflux::PenaltyCreatePVector(Vec pg,double kmax){
     384void  Numericalflux::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){
    385385
    386386        /*No penalty loads applied, do nothing: */
    … …  
    417417
    418418        GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES_INTERNAL);
    419         GetDofList(&doflist);
     419        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    420420
    421421        /*Retrieve all inputs and parameters we will be needing: */
    … …  
    523523        }
    524524
    525         GetDofList(&doflist);
     525        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    526526
    527527        /* Start  looping on the number of gaussian points: */
    … …  
    617617        }
    618618
    619         GetDofList(&doflist);
     619        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    620620
    621621        /* Start  looping on the number of gaussian points: */
    … …  
    646646/*}}}*/
    647647/*FUNCTION Numericalflux::GetDofList {{{1*/
    648 void  Numericalflux::GetDofList(int** pdoflist){
     648void  Numericalflux::GetDofList(int** pdoflist,int approximation,int setenum){
    649649
    650650        int i,j;
    … …  
    671671        /*Figure out size of doflist: */
    672672        for(i=0;i<numberofnodes;i++){
    673                 numberofdofs+=nodes[i]->GetNumberOfDofs();
     673                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    674674        }
    675675
    … …  
    680680        count=0;
    681681        for(i=0;i<numberofnodes;i++){
    682                 nodes[i]->GetDofList(doflist+count);
    683                 count+=nodes[i]->GetNumberOfDofs();
     682                nodes[i]->GetDofList(doflist+count,approximation,setenum);
     683                count+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    684684        }
    685685

  • issm/trunk/src/c/objects/Loads/Numericalflux.h

    r5739 r5772  
    6565                void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    6666                void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    67                 void  CreateKMatrix(Mat Kgg);
    68                 void  CreatePVector(Vec pg);
    69                 void  PenaltyCreateKMatrix(Mat Kgg,double kmax);
    70                 void  PenaltyCreatePVector(Vec pg,double kmax);
     67                void  CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     68                void  CreatePVector(Vec pg, Vec pf);
     69                void  PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax);
     70                void  PenaltyCreatePVector(Vec pg,Vec pf, double kmax);
    7171                bool  InAnalysis(int analysis_type);
    7272                /*}}}*/
    7373                /*Numericalflux management:{{{1*/
    74                 void  GetDofList(int** pdoflist);
     74                void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
    7575                void  GetNormal(double* normal,double xyz_list[4][3]);
    7676                void  CreateKMatrixInternal(Mat Kgg);

  • issm/trunk/src/c/objects/Loads/Pengrid.cpp

    r5735 r5772  
    281281/*}}}1*/
    282282/*FUNCTION Pengrid::CreateKMatrix {{{1*/
    283 void  Pengrid::CreateKMatrix(Mat Kgg){
     283void  Pengrid::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){
    284284
    285285        /*No loads applied, do nothing: */
    … …  
    289289/*}}}1*/
    290290/*FUNCTION Pengrid::CreatePVector {{{1*/
    291 void  Pengrid::CreatePVector(Vec pg){
     291void  Pengrid::CreatePVector(Vec pg,Vec pf){
    292292
    293293        /*No loads applied, do nothing: */
    … …  
    297297/*}}}1*/
    298298/*FUNCTION Pengrid::PenaltyCreateMatrix {{{1*/
    299 void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,double kmax){
     299void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){
    300300
    301301        int analysis_type;
    … …  
    320320/*}}}1*/
    321321/*FUNCTION Pengrid::PenaltyCreatePVector {{{1*/
    322 void  Pengrid::PenaltyCreatePVector(Vec pg,double kmax){
     322void  Pengrid::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){
    323323
    324324        int analysis_type;
    … …  
    419419/*Pengrid management:*/
    420420/*FUNCTION Pengrid::GetDofList {{{1*/
    421 void  Pengrid::GetDofList(int** pdoflist){
     421void  Pengrid::GetDofList(int** pdoflist,int approximation,int setenum){
    422422
    423423        int  i,j;
    … …  
    431431
    432432        /*First, figure out size of doflist: */
    433         numberofdofs=node->GetNumberOfDofs();
     433        numberofdofs=node->GetNumberOfDofs(approximation,setenum);
    434434
    435435        /*Allocate: */
    … …  
    437437
    438438        /*Populate: */
    439         node->GetDofList(doflist);
     439        node->GetDofList(doflist,approximation,setenum);
    440440
    441441        /*Assign output pointers:*/
    … …  
    567567
    568568        /*Get dof list: */
    569         GetDofList(&doflist);
     569        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    570570
    571571        //recover slope: */
    … …  
    620620
    621621        /*Get dof list: */
    622         GetDofList(&doflist);
     622        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    623623       
    624624        //Compute pressure melting point
    … …  
    656656
    657657        /*Get dof list: */
    658         GetDofList(&doflist);
     658        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    659659
    660660        Ke[0][0]=kmax*pow((double)10,penalty_offset);
    … …  
    693693
    694694        /*Get dof list: */
    695         GetDofList(&doflist);
     695        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    696696
    697697        //First recover pressure and temperature values, using the element: */
    … …  
    756756
    757757        /*Get dof list: */
    758         GetDofList(&doflist);
     758        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    759759
    760760        //First recover pressure  and penalty_offset

  • issm/trunk/src/c/objects/Loads/Pengrid.h

    r5730 r5772  
    7171                void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    7272                void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    73                 void  CreateKMatrix(Mat Kgg);
    74                 void  CreatePVector(Vec pg);
    75                 void  PenaltyCreateKMatrix(Mat Kgg,double kmax);
    76                 void  PenaltyCreatePVector(Vec pg,double kmax);
     73                void  CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     74                void  CreatePVector(Vec pg, Vec pf);
     75                void  PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax);
     76                void  PenaltyCreatePVector(Vec pg,Vec pf, double kmax);
    7777                bool  InAnalysis(int analysis_type);
    7878                /*}}}*/
    7979                /*Pengrid management {{{1*/
    8080                void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax);
    81                 void  GetDofList(int** pdoflist);
     81                void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
    8282                void  PenaltyCreateKMatrixThermal(Mat Kgg,double kmax);
    8383                void  PenaltyCreateKMatrixMelting(Mat Kgg,double kmax);

  • issm/trunk/src/c/objects/Loads/Penpair.cpp

    r5596 r5772  
    184184/*}}}1*/
    185185/*FUNCTION Penpair::CreateKMatrix {{{1*/
    186 void  Penpair::CreateKMatrix(Mat Kgg){
     186void  Penpair::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){
    187187        /*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/
    188188        /*No loads applied, do nothing: */
    … …  
    192192/*}}}1*/
    193193/*FUNCTION Penpair::CreatePVector {{{1*/
    194 void  Penpair::CreatePVector(Vec pg){
     194void  Penpair::CreatePVector(Vec pg,Vec pf){
    195195
    196196        /*No loads applied, do nothing: */
    … …  
    200200/*}}}1*/
    201201/*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/
    202 void  Penpair::PenaltyCreateKMatrix(Mat Kgg,double kmax){
     202void  Penpair::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){
    203203        int analysis_type;
    204204
    … …  
    251251/*}}}1*/
    252252/*FUNCTION Penpair::PenaltyCreatePVector {{{1*/
    253 void  Penpair::PenaltyCreatePVector(Vec pg,double kmax){
     253void  Penpair::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){
    254254        /*No loads applied, do nothing: */
    255255        return;
    … …  
    297297
    298298        /*Get dof list: */
    299         GetDofList(&doflist);
     299        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    300300
    301301        /*recover pointers: */
    … …  
    340340
    341341        /*Get dof list: */
    342         GetDofList(&doflist);
     342        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    343343
    344344        /*recover pointers: */
    … …  
    378378/*}}}1*/
    379379/*FUNCTION Penpair::GetDofList {{{1*/
    380 void  Penpair::GetDofList(int** pdoflist){
     380void  Penpair::GetDofList(int** pdoflist,int approximation,int setenum){
    381381
    382382        int i,j;
    … …  
    398398        /*First, figure out size of doflist: */
    399399        for(i=0;i<2;i++){
    400                 numberofdofs+=nodes[i]->GetNumberOfDofs();
     400                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    401401        }
    402402
    … …  
    407407        count=0;
    408408        for(i=0;i<2;i++){
    409                 nodes[i]->GetDofList(doflist+count);
    410                 count+=nodes[i]->GetNumberOfDofs();
     409                nodes[i]->GetDofList(doflist+count,approximation,setenum);
     410                count+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    411411        }
    412412

  • issm/trunk/src/c/objects/Loads/Penpair.h

    r5510 r5772  
    5757                void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    5858                void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    59                 void  CreateKMatrix(Mat Kgg);
    60                 void  CreatePVector(Vec pg);
    61                 void  PenaltyCreateKMatrix(Mat Kgg,double kmax);
    62                 void  PenaltyCreatePVector(Vec pg,double kmax);
     59                void  CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     60                void  CreatePVector(Vec pg, Vec pf);
     61                void  PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax);
     62                void  PenaltyCreatePVector(Vec pg,Vec pf, double kmax);
    6363                bool  InAnalysis(int analysis_type);
    6464                /*}}}*/
    … …  
    6666                void  PenaltyCreateKMatrixDiagnosticHoriz(Mat Kgg,double kmax);
    6767                void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax);
    68                 void  GetDofList(int** pdoflist);
     68                void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
    6969                /*}}}*/
    7070};

  • issm/trunk/src/c/objects/Loads/Riftfront.cpp

    r5655 r5772  
    376376/*}}}*/
    377377/*FUNCTION Riftfront::CreateKMatrix {{{1*/
    378 void  Riftfront::CreateKMatrix(Mat Kgg){
     378void  Riftfront::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){
    379379        /*do nothing: */
    380380}
    381381/*}}}1*/
    382382/*FUNCTION Riftfront::CreatePVector {{{1*/
    383 void  Riftfront::CreatePVector(Vec pg){
     383void  Riftfront::CreatePVector(Vec pg,Vec pf){
    384384        /*do nothing: */
    385385}
    386386/*}}}1*/
    387387/*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/
    388 void  Riftfront::PenaltyCreateKMatrix(Mat Kgg,double kmax){
     388void  Riftfront::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){
    389389
    390390        int         i;
    … …  
    419419       
    420420        /* Get node coordinates and dof list: */
    421         GetDofList(&doflist);
     421        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    422422
    423423        /* Set Ke_gg to 0: */
    … …  
    499499/*}}}1*/
    500500/*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/
    501 void  Riftfront::PenaltyCreatePVector(Vec pg,double kmax){
     501void  Riftfront::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){
    502502
    503503        int         i                     ,j;
    … …  
    544544
    545545        /* Get node coordinates and dof list: */
    546         GetDofList(&doflist);
     546        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    547547
    548548        /*Get some inputs: */
    … …  
    731731/*}}}1*/
    732732/*FUNCTION Riftfront::GetDofList {{{1*/
    733 void  Riftfront::GetDofList(int** pdoflist){
     733void  Riftfront::GetDofList(int** pdoflist,int approximation,int setenum){
    734734
    735735        int i,j;
    … …  
    751751        /*Figure out size of doflist: */
    752752        for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
    753                 numberofdofs+=nodes[i]->GetNumberOfDofs();
     753                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    754754        }
    755755
    … …  
    760760        count=0;
    761761        for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
    762                 nodes[i]->GetDofList(doflist+count);
    763                 count+=nodes[i]->GetNumberOfDofs();
     762                nodes[i]->GetDofList(doflist+count,approximation,setenum);
     763                count+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    764764        }
    765765

  • issm/trunk/src/c/objects/Loads/Riftfront.h

    r5655 r5772  
    7474                void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    7575                void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    76                 void  CreateKMatrix(Mat Kgg);
    77                 void  CreatePVector(Vec pg);
    78                 void  PenaltyCreateKMatrix(Mat Kgg,double kmax);
    79                 void  PenaltyCreatePVector(Vec pg,double kmax);
     76                void  CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
     77                void  CreatePVector(Vec pg, Vec pf);
     78                void  PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax);
     79                void  PenaltyCreatePVector(Vec pg,Vec pf, double kmax);
    8080                bool  InAnalysis(int analysis_type);
    8181                /*}}}*/
    8282                /*Riftfront specific routines: {{{1*/
    83                 void  GetDofList(int** pdoflist);
     83                void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
    8484                bool  PreStable();
    8585                void  SetPreStable();

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

    r5578 r5772  
    6363        /*Intermediary*/
    6464        int k;
    65         int numdofs;
     65        int gsize;
    6666
    6767        /*id: */
    … …  
    7272        /*indexing:*/
    7373        DistributeNumDofs(&this->indexing,analysis_type,iomodel->vertices_type+io_index); //number of dofs per node
    74         numdofs=this->indexing.numberofdofs;
     74        gsize=this->indexing.gsize;
    7575
    7676        /*Hooks*/
    … …  
    9595                        if (!iomodel->vertices_type) ISSMERROR("iomodel->vertices_type is NULL");
    9696                        if (iomodel->vertices_type[io_index]==MacAyealApproximationEnum && !iomodel->gridonbed[io_index]){
    97                                 for(k=1;k<=numdofs;k++) this->FreezeDof(k);
     97                                for(k=1;k<=gsize;k++) this->FreezeDof(k);
    9898                        }
    9999                        if (iomodel->vertices_type[io_index]==MacAyealPattynApproximationEnum && iomodel->gridonmacayeal[io_index]){
    100100                                if(!iomodel->gridonbed[io_index]){
    101                                         for(k=1;k<=numdofs;k++) this->FreezeDof(k);
     101                                        for(k=1;k<=gsize;k++) this->FreezeDof(k);
    102102                                }
    103103                        }
    … …  
    106106                if (!iomodel->gridonhutter) ISSMERROR("iomodel->gridonhutter is NULL");
    107107                if (iomodel->gridonhutter[io_index]){
    108                         for(k=1;k<=numdofs;k++){
     108                        for(k=1;k<=gsize;k++){
    109109                                this->FreezeDof(k);
    110110                        }
    … …  
    117117                if (!iomodel->gridonhutter) ISSMERROR("iomodel->gridonhutter is NULL");
    118118                if (!iomodel->gridonhutter[io_index]){
    119                         for(k=1;k<=numdofs;k++){
     119                        for(k=1;k<=gsize;k++){
    120120                                this->FreezeDof(k);
    121121                        }
    … …  
    136136                        if (!iomodel->gridonbed) ISSMERROR("iomodel->gridonbed is NULL");
    137137                        if (!iomodel->gridonbed[io_index]){
    138                                 for(k=1;k<=numdofs;k++){
     138                                for(k=1;k<=gsize;k++){
    139139                                        this->FreezeDof(k);
    140140                                }
    … …  
    199199        int   enum_type=0;
    200200        char* marshalled_inputs=NULL;
    201         int   marshalled_inputs_size;
     201        int   marshalled_inputssize;
    202202
    203203        /*recover marshalled_dataset: */
    … …  
    220220
    221221        /*Marshall inputs: */
    222         marshalled_inputs_size=inputs->MarshallSize();
     222        marshalled_inputssize=inputs->MarshallSize();
    223223        marshalled_inputs=inputs->Marshall();
    224         memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
    225         marshalled_dataset+=marshalled_inputs_size;
     224        memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputssize*sizeof(char));
     225        marshalled_dataset+=marshalled_inputssize;
    226226
    227227        /*Free ressources:*/
    … …  
    294294
    295295        /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
    296          * datasets, using internal ids and offsets hidden in hooks: */
     296         * datasets, using internal ids and off_sets hidden in hooks: */
    297297        hvertex->configure(verticesin);
    298298
    … …  
    303303}
    304304/*FUNCTION Node::GetDof {{{1*/
    305 int   Node::GetDof(int dofindex){
    306 
    307         return indexing.doflist[dofindex];
     305int   Node::GetDof(int dofindex,int setenum){
     306
     307        if(setenum==GsetEnum){
     308                return indexing.gdoflist[dofindex];
     309        }
     310        else if(setenum==FsetEnum){
     311                return indexing.fdoflist[dofindex];
     312        }
     313        else if(setenum==SsetEnum){
     314                return indexing.sdoflist[dofindex];
     315        }
     316        else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
    308317
    309318}
    … …  
    320329/*}}}*/
    321330/*FUNCTION Node::GetDofList{{{1*/
    322 void  Node::GetDofList(int* outdoflist,int approximation_enum){
     331void  Node::GetDofList(int* outdoflist,int approximation_enum,int setenum){
    323332        int i;
    324333        int count=0;
    325 
    326         if(approximation_enum && this->indexing.doftype){ //if an approximation is specified and this node has several approsimations
    327                 for(i=0;i<this->indexing.numberofdofs;i++){
    328                         if(indexing.doftype[i]==approximation_enum){
    329                                 outdoflist[count]=indexing.doflist[i];
    330                                 count++;
    331                         }
    332                 }
    333                 ISSMASSERT(count);
     334               
     335        if(approximation_enum==NoneApproximationEnum){
     336                if(setenum==GsetEnum)for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i];
     337                if(setenum==FsetEnum)for(i=0;i<this->indexing.fsize;i++) outdoflist[i]=indexing.fdoflist[i];
     338                if(setenum==SsetEnum)for(i=0;i<this->indexing.ssize;i++) outdoflist[i]=indexing.sdoflist[i];
    334339        }
    335340        else{
    336                 for(i=0;i<this->indexing.numberofdofs;i++){
    337                         outdoflist[i]=indexing.doflist[i];
    338                 }
     341
     342                if(setenum==GsetEnum){
     343                        if(indexing.doftype){
     344                                count=0;
     345                                for(i=0;i<this->indexing.gsize;i++){
     346                                        if(indexing.doftype[i]==approximation_enum){
     347                                                outdoflist[count]=indexing.gdoflist[i];
     348                                                count++;
     349                                        }
     350                                }
     351                                ISSMASSERT(count);
     352                        }
     353                        else for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i];
     354                }
     355                else if(setenum==FsetEnum){
     356                        if(indexing.doftype){
     357                                count=0;
     358                                for(i=0;i<this->indexing.gsize;i++){
     359                                        if((indexing.doftype[i]==approximation_enum) && (indexing.f_set[i])){
     360                                                outdoflist[count]=indexing.fdoflist[i];
     361                                                count++;
     362                                        }
     363                                }
     364                                ISSMASSERT(count);
     365                        }
     366                        else for(i=0;i<this->indexing.fsize;i++) outdoflist[i]=indexing.fdoflist[i];
     367                }
     368                else if(setenum==SsetEnum){
     369                        if(indexing.doftype){
     370                                count=0;
     371                                for(i=0;i<this->indexing.gsize;i++){
     372                                        if((indexing.doftype[i]==approximation_enum) && (indexing.s_set[i])){
     373                                                outdoflist[count]=indexing.sdoflist[i];
     374                                                count++;
     375                                        }
     376                                }
     377                                ISSMASSERT(count);
     378                        }
     379                        else for(i=0;i<this->indexing.ssize;i++) outdoflist[i]=indexing.sdoflist[i];
     380                }
     381                else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
     382        }
     383}
     384/*}}}*/
     385/*FUNCTION Node::GetInternalDofList{{{1*/
     386void  Node::GetInternalDofList(int* outdoflist,int approximation_enum,int setenum){
     387        int i;
     388        int count=0;
     389        int count2=0;
     390               
     391        if(approximation_enum==NoneApproximationEnum){
     392                if(setenum==GsetEnum)for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=i;
     393                else if(setenum==FsetEnum){
     394                        count=0;
     395                        for(i=0;i<this->indexing.fsize;i++){
     396                                if(indexing.f_set[i]){
     397                                        outdoflist[count]=i;
     398                                        count++;
     399                                }
     400                        }
     401                }
     402                else if(setenum==SsetEnum){
     403                        count=0;
     404                        for(i=0;i<this->indexing.ssize;i++){
     405                                if(indexing.s_set[i]){
     406                                        outdoflist[count]=i;
     407                                        count++;
     408                                }
     409                        }
     410                }
     411                else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
     412        }
     413        else{
     414
     415                if(setenum==GsetEnum){
     416                        if(indexing.doftype){
     417                                count=0;
     418                                for(i=0;i<this->indexing.gsize;i++){
     419                                        if(indexing.doftype[i]==approximation_enum){
     420                                                outdoflist[count]=count;
     421                                                count++;
     422                                        }
     423                                }
     424                                ISSMASSERT(count);
     425                        }
     426                        else for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=i;
     427                }
     428                else if(setenum==FsetEnum){
     429
     430                        if(indexing.doftype){
     431                                count=0;
     432                                count2=0;
     433                                for(i=0;i<this->indexing.gsize;i++){
     434                                        if((indexing.doftype[i]==approximation_enum) && (indexing.f_set[i])){
     435                                                outdoflist[count]=count2;
     436                                                count++;
     437                                        }
     438                                        if(indexing.doftype[i]==approximation_enum)count2++;
     439                                }
     440                                ISSMASSERT(count);
     441                        }
     442                        else{
     443
     444                                count=0;
     445                                for(i=0;i<this->indexing.gsize;i++){
     446                                        if(indexing.f_set[i]){
     447                                                outdoflist[count]=i;
     448                                                count++;
     449                                        }
     450                                }
     451                        }
     452                }
     453                else if(setenum==SsetEnum){
     454                        if(indexing.doftype){
     455                                count=0;
     456                                count2=0;
     457                                for(i=0;i<this->indexing.gsize;i++){
     458                                        if((indexing.doftype[i]==approximation_enum) && (indexing.s_set[i])){
     459                                                outdoflist[count]=count2;
     460                                                count++;
     461                                        }
     462                                        if(indexing.doftype[i]==approximation_enum)count2++;
     463                                }
     464                                ISSMASSERT(count);
     465                        }
     466                        else{
     467                                count=0;
     468                                for(i=0;i<this->indexing.gsize;i++){
     469                                        if(indexing.s_set[i]){
     470                                                outdoflist[count]=i;
     471                                                count++;
     472                                        }
     473                                }
     474                        }
     475                }
     476                else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
    339477        }
    340478}
    … …  
    380518/*Node numerics:*/
    381519/*FUNCTION Node::ApplyConstraints{{{1*/
    382 void  Node::ApplyConstraint(Vec yg,int dof,double value){
     520void  Node::ApplyConstraint(int dof,double value){
    383521
    384522        int index;
    385523
    386         /*First, dof should be added in the s set, describing which
     524        /*Dof should be added in the s set, describing which
    387525         * dofs are constrained to a certain value (dirichlet boundary condition*/
    388 
    389526        DofInSSet(dof-1);
    390 
    391         /*Second, we should add value into yg, at dof corresponding to doflist[dof], unless
    392          *  we are a clone!*/
    393 
    394         if(!indexing.clone){
    395 
    396                 index=indexing.doflist[dof-1]; //matlab indexing
    397 
    398                 VecSetValues(yg,1,&index,&value,INSERT_VALUES);
    399 
    400         }
    401 
     527        this->indexing.svalues[dof-1]=value;
    402528}
    403529/*}}}*/
    … …  
    410536        int i;
    411537
    412         for(i=0;i<this->indexing.numberofdofs;i++){
     538        for(i=0;i<this->indexing.gsize;i++){
    413539
    414540                /*g set: */
    415                 VecSetValues(pv_g,1,&indexing.doflist[i],&gvalue,INSERT_VALUES);
     541                VecSetValues(pv_g,1,&indexing.gdoflist[i],&gvalue,INSERT_VALUES);
    416542               
    417543                /*f set: */
    418544                value=(double)this->indexing.f_set[i];
    419                 VecSetValues(pv_f,1,&indexing.doflist[i],&value,INSERT_VALUES);
     545                VecSetValues(pv_f,1,&indexing.gdoflist[i],&value,INSERT_VALUES);
    420546
    421547                /*s set: */
    422548                value=(double)this->indexing.s_set[i];
    423                 VecSetValues(pv_s,1,&indexing.doflist[i],&value,INSERT_VALUES);
    424 
    425         }
     549                VecSetValues(pv_s,1,&indexing.gdoflist[i],&value,INSERT_VALUES);
     550
     551        }
     552
     553
     554}
     555/*}}}*/
     556/*FUNCTION Node::CreateNodalConstraints{{{1*/
     557void  Node::CreateNodalConstraints(Vec ys){
     558
     559        int i;
     560        double* values=NULL;
     561        int count;
     562
     563        /*Recover values for s set and plug them in constraints vector: */
     564        if(this->indexing.ssize){
     565                values=(double*)xmalloc(this->indexing.ssize*sizeof(double));
     566                count=0;
     567                for(i=0;i<this->indexing.gsize;i++){
     568                        if(this->indexing.s_set[i])values[count]=this->indexing.svalues[i];
     569                        count++;
     570                }
     571
     572                /*Add values into constraint vector: */
     573                VecSetValues(ys,this->indexing.ssize,this->indexing.sdoflist,values,INSERT_VALUES);
     574        }
     575
     576        /*Free ressources:*/
     577        xfree((void**)&values);
    426578
    427579
    … …  
    467619/*}}}*/
    468620/*FUNCTION Node::GetNumberOfDofs{{{1*/
    469 int   Node::GetNumberOfDofs(int approximation_enum){
     621int   Node::GetNumberOfDofs(int approximation_enum,int setenum){
     622
     623        /*Get number of degrees of freedom in a node, for a certain set (g,f or s-set)
     624         *and for a certain approximation type: */
    470625       
    471626        int i;
    472627        int numdofs=0;
    473628
    474         /*Count the dofs if an approximation is specified and the node contains several dofs type*/
    475         if(approximation_enum && this->indexing.doftype){
    476                 for(i=0;i<this->indexing.numberofdofs;i++){
    477                         if(this->indexing.doftype[i]==approximation_enum) numdofs++;
    478                 }
    479         }
    480         else numdofs=this->indexing.numberofdofs;
    481 
    482         ISSMASSERT(numdofs);
     629        if(approximation_enum==NoneApproximationEnum){
     630                if (setenum==GsetEnum) numdofs=this->indexing.gsize;
     631                else if (setenum==FsetEnum) numdofs=this->indexing.fsize;
     632                else if (setenum==SsetEnum) numdofs=this->indexing.ssize;
     633                else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
     634        }
     635        else{
     636                if(setenum==GsetEnum){
     637                        if(this->indexing.doftype){
     638                                numdofs=0;
     639                                for(i=0;i<this->indexing.gsize;i++){
     640                                        if(this->indexing.doftype[i]==approximation_enum) numdofs++;
     641                                }
     642                        }
     643                        else numdofs=this->indexing.gsize;
     644                }
     645                else if (setenum==FsetEnum){
     646                        if(this->indexing.doftype){
     647                                numdofs=0;
     648                                for(i=0;i<this->indexing.gsize;i++){
     649                                        if((this->indexing.doftype[i]==approximation_enum) && (this->indexing.f_set[i])) numdofs++;
     650                                }
     651                        }
     652                        else numdofs=this->indexing.fsize;
     653                }
     654                else if (setenum==SsetEnum){
     655                        if(this->indexing.doftype){
     656                                numdofs=0;
     657                                for(i=0;i<this->indexing.gsize;i++){
     658                                        if((this->indexing.doftype[i]==approximation_enum) && (this->indexing.s_set[i])) numdofs++;
     659                                }
     660                        }
     661                        else numdofs=this->indexing.ssize;
     662                }
     663                else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
     664        }
    483665        return numdofs;
    484 
    485666}
    486667/*}}}*/
    … …  
    625806/* DofObject routines:*/
    626807/*FUNCTION Node::DistributeDofs{{{1*/
    627 void  Node::DistributeDofs(int* pdofcount){
     808void  Node::DistributeDofs(int* pdofcount,int setenum){
    628809
    629810        int i;
    … …  
    638819        }
    639820
    640         /*This node should distribute dofs, go ahead: */
    641         for(i=0;i<this->indexing.numberofdofs;i++){
    642                 indexing.doflist[i]=dofcount+i;
    643         }
    644         dofcount+=this->indexing.numberofdofs;
     821        /*This node should distribute dofs for setenum set (eg, f_set or s_set), go ahead: */
     822
     823        if(setenum==GsetEnum){
     824                for(i=0;i<this->indexing.gsize;i++){
     825                        indexing.gdoflist[i]=dofcount+i;
     826                }
     827                dofcount+=this->indexing.gsize;
     828        }
     829        else if(setenum==FsetEnum){
     830                this->indexing.InitSet(setenum);
     831                for(i=0;i<this->indexing.fsize;i++){
     832                        indexing.fdoflist[i]=dofcount+i;
     833                }
     834                dofcount+=this->indexing.fsize;
     835        }
     836        else if(setenum==SsetEnum){
     837                this->indexing.InitSet(setenum);
     838                for(i=0;i<this->indexing.ssize;i++){
     839                        indexing.sdoflist[i]=dofcount+i;
     840                }
     841                dofcount+=this->indexing.ssize;
     842        }
     843        else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
     844
    645845
    646846        /*Assign output pointers: */
    … …  
    649849}
    650850/*}}}*/
    651 /*FUNCTION Node::OffsetDofs{{{1*/
    652 void  Node::OffsetDofs(int dofcount){
     851/*FUNCTION Node::Off_setDofs{{{1*/
     852void  Node::OffsetDofs(int dofcount,int setenum){
    653853       
    654854        int i;
    … …  
    656856       
    657857        if(indexing.clone){
    658                 /*This node is a clone, don't offset the dofs!: */
     858                /*This node is a clone, don't off_set the dofs!: */
    659859                return;
    660860        }
    661861
    662         /*This node should offset the dofs, go ahead: */
    663         for(i=0;i<this->indexing.numberofdofs;i++){
    664                 indexing.doflist[i]+=dofcount;
    665         }
     862        /*This node should off_set the dofs, go ahead: */
     863        if(setenum==GsetEnum){
     864                for(i=0;i<this->indexing.gsize;i++) indexing.gdoflist[i]+=dofcount;
     865        }
     866        else if(setenum==FsetEnum){
     867                for(i=0;i<this->indexing.fsize;i++) indexing.fdoflist[i]+=dofcount;
     868        }
     869        else if(setenum==SsetEnum){
     870                for(i=0;i<this->indexing.ssize;i++) indexing.sdoflist[i]+=dofcount;
     871        }
     872        else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
    666873}
    667874/*}}}*/
    668875/*FUNCTION Node::ShowTrueDofs{{{1*/
    669 void  Node::ShowTrueDofs(int* truedofs, int ncols){
     876void  Node::ShowTrueDofs(int* truedofs, int ncols,int setenum){
    670877
    671878        int j;
    … …  
    676883
    677884        /*Ok, we are not a clone, just plug our dofs into truedofs: */
    678         for(j=0;j<this->indexing.numberofdofs;j++){
    679                 *(truedofs+ncols*sid+j)=indexing.doflist[j];
    680         }
     885        if(setenum==GsetEnum)for(j=0;j<this->indexing.gsize;j++)  *(truedofs+ncols*sid+j)=indexing.gdoflist[j];
     886        else if(setenum==FsetEnum)for(j=0;j<this->indexing.fsize;j++)  *(truedofs+ncols*sid+j)=indexing.fdoflist[j];
     887        else if(setenum==SsetEnum)for(j=0;j<this->indexing.ssize;j++)  *(truedofs+ncols*sid+j)=indexing.sdoflist[j];
     888        else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
     889
    681890}
    682891/*}}}*/
    683892/*FUNCTION Node::UpdateCloneDofs{{{1*/
    684 void  Node::UpdateCloneDofs(int* alltruedofs,int ncols){
     893void  Node::UpdateCloneDofs(int* alltruedofs,int ncols,int setenum){
    685894
    686895        int j;
    … …  
    693902        /*Ok, we are a clone node, but we did not create the dofs for this node.
    694903         *      * Therefore, our doflist is garbage right now. Go pick it up in the alltruedofs: */
    695         for(j=0;j<this->indexing.numberofdofs;j++){
    696                 indexing.doflist[j]=*(alltruedofs+ncols*sid+j);
    697         }
     904        if(setenum==GsetEnum)for(j=0;j<this->indexing.gsize;j++) indexing.gdoflist[j]=*(alltruedofs+ncols*sid+j);
     905        else if(setenum==FsetEnum)for(j=0;j<this->indexing.fsize;j++) indexing.fdoflist[j]=*(alltruedofs+ncols*sid+j);
     906        else if(setenum==SsetEnum)for(j=0;j<this->indexing.ssize;j++) indexing.sdoflist[j]=*(alltruedofs+ncols*sid+j);
     907        else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!");
    698908
    699909}

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

    r5557 r5772  
    6565                /*Node numerical routines {{{1*/
    6666                void  Configure(DataSet* nodes,Vertices* vertices);
     67                void  CreateNodalConstraints(Vec ys);
    6768                void  SetCurrentConfiguration(DataSet* nodes,Vertices* vertices);
    6869                int   Sid(void);
    … …  
    7273                bool  InAnalysis(int analysis_type);
    7374                int   GetApproximation();
    74                 int   GetNumberOfDofs(int approximation_enum=0);
     75                int   GetNumberOfDofs(int approximation_enum,int setenum);
    7576                int   IsClone();
    76                 void  ApplyConstraint(Vec yg,int dof,double value);
     77                void  ApplyConstraint(int dof,double value);
    7778                void  DofInSSet(int dof);
    78                 int   GetDof(int dofindex);
     79                int   GetDof(int dofindex,int setenum);
    7980                void  CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s);
    8081                int   GetConnectivity();
    81                 void  GetDofList(int* poutdoflist,int approximation_enum=0);
     82                void  GetDofList(int* poutdoflist,int approximation_enum,int setenum);
     83                void  GetInternalDofList(int* poutdoflist,int approximation_enum,int setenum);
    8284                int   GetDofList1(void);
    8385                double GetX();
    … …  
    9294                /*}}}*/
    9395                /*Dof Object routines {{{1*/
    94                 void  DistributeDofs(int* pdofcount);
    95                 void  OffsetDofs(int dofcount);
    96                 void  ShowTrueDofs(int* truerows,int ncols);
    97                 void  UpdateCloneDofs(int* alltruerows,int ncols);
     96                void  DistributeDofs(int* pdofcount,int setenum);
     97                void  OffsetDofs(int dofcount,int setenum);
     98                void  ShowTrueDofs(int* truerows,int ncols,int setenum);
     99                void  UpdateCloneDofs(int* alltruerows,int ncols,int setenum);
    98100                void  SetClone(int* minranks);
    99101                void  CreatePartition(Vec partition);

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

    r5740 r5772  
    7373#include "./Materials/Matpar.h"
    7474
     75/*Numerics:*/
     76#include "./Numerics/ElementMatrix.h"
     77#include "./Numerics/ElementVector.h"
     78
    7579
    7680/*Params: */

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

    r4445 r5772  
    1010       
    1111        int verbose=0;
    12         Vec ug=NULL;
     12        Vec yg=NULL;
     13        Vec ys=NULL;
    1314        int analysis_counter;
    1415                       
    … …  
    1920        femmodel->SetCurrentConfiguration(analysis_type);
    2021
    21         GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,  femmodel->parameters);
     22        GetSolutionFromInputsx( &yg, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,  femmodel->parameters);
    2223
    23         /*For this analysis_type, free existing boundary condition vectors: */
     24        /*For this analysis_type, free existing boundary condition vector: */
    2425        analysis_counter=femmodel->analysis_counter;
    25 
    26         //global dof set
    27         VecFree(&femmodel->m_yg[analysis_counter]);
    28         //in the s-set
    2926        VecFree(&femmodel->m_ys[analysis_counter]);
    3027
    31         //Now, duplicate ug (the solution vector) into the boundary conditions vector on the g-set
    32         VecDuplicatePatch(&femmodel->m_yg[analysis_counter],ug);
    33        
    3428        //Reduce from g to s set
    35         Reducevectorgtosx(&femmodel->m_ys[analysis_counter],femmodel->m_yg[analysis_counter],femmodel->m_nodesets[analysis_counter]);
     29        Reducevectorgtosx(&ys,yg,femmodel->m_nodesets[analysis_counter]);
     30
     31        /*Plug into femmodel->m_ys: */
     32        femmodel->m_ys[analysis_counter]=ys;
    3633
    3734        /*Free ressources:*/
    38         VecFree(&ug);
    39 
     35        VecFree(&yg);
    4036}

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

    r5695 r5772  
    1717        Vec pg  = NULL, pf=NULL;
    1818
    19         SystemMatricesx(&Kgg, &pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     19        SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    2020
    2121        Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);

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

    r5695 r5772  
    1212void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads){
    1313
    14         /*intermediary: */
    15         Mat Kgg = NULL, Kff = NULL, Kfs   = NULL;
    16         Vec ug  = NULL, uf  = NULL, old_ug= NULL, old_uf = NULL;
    17         Vec pg  = NULL, pf  = NULL;
    18         Loads* loads=NULL;
    19         int converged;
    20         int constraints_converged;
    21         int num_unstable_constraints;
    22         int count;
    23         int min_mechanical_constraints;
    24         int max_nonlinear_iterations;
    25 
    2614        /*parameters:*/
    27         int verbose=0;
     15        bool kff=false;
    2816
    2917        /*Recover parameters: */
    30         femmodel->parameters->FindParam(&verbose,VerboseEnum);
    31         femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum);
    32         femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum);
    33        
    34         /*Were loads requested as output? : */
    35         if(conserve_loads){
    36                 loads=(Loads*)femmodel->loads->Copy(); //protect loads from being modified by the solution
    37         }
    38         else{
    39                 loads=(Loads*)femmodel->loads; //modify loads  in this solution
    40         }
     18        femmodel->parameters->FindParam(&kff,KffEnum);
    4119
    42         count=1;
    43         converged=0;
     20        if(!kff) solver_diagnostic_nonlinear_kgg(femmodel,conserve_loads);
     21        else     solver_diagnostic_nonlinear_kff(femmodel,conserve_loads);
    4422
    45         /*Start non-linear iteration using input velocity: */
    46         GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters);
    47         Reducevectorgtofx(&uf, ug, femmodel->nodesets);
    48 
    49         //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
    50         InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
    51 
    52         for(;;){
    53 
    54                 //save pointer to old velocity
    55                 VecFree(&old_ug);old_ug=ug;
    56                 VecFree(&old_uf);old_uf=uf;
    57 
    58                 SystemMatricesx(&Kgg, &pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
    59                
    60                 Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);
    61        
    62                 Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters); VecFree(&pg); MatFree(&Kfs);
    63 
    64                 Solverx(&uf, Kff, pf, old_uf, femmodel->parameters);
    65 
    66                 Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters);
    67 
    68                 InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
    69 
    70                 PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
    71                 if(verbose)_printf_("   number of unstable constraints: %i\n",num_unstable_constraints);
    72 
    73                 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf);
    74                
    75                 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
    76 
    77                 //rift convergence
    78                 if (!constraints_converged) {
    79                         if (converged){
    80                                 if (num_unstable_constraints <= min_mechanical_constraints) converged=1;
    81                                 else converged=0;
    82                         }
    83                 }
    84 
    85                 /*Increase count: */
    86                 count++;
    87                 if(converged==1)break;
    88                 if(count>=max_nonlinear_iterations){
    89                         _printf_("   maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
    90                         break;
    91                 }
    92         }
    93 
    94         /*clean-up*/
    95         if(conserve_loads) delete loads;
    96         VecFree(&uf);
    97         VecFree(&ug);
    98         VecFree(&old_uf);
    99         VecFree(&old_ug);
    100        
    10123}

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

    r5695 r5772  
    1515        Vec pg  = NULL, pf  = NULL;
    1616
    17         SystemMatricesx(&Kgg,&pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     17        SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    1818
    1919        Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);

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

    r5707 r5772  
    5151
    5252
    53                 SystemMatricesx(&Kgg,&pg,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     53                SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    5454
    5555                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);

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

    r4839 r5772  
    1414void solver_thermal_nonlinear(FemModel* femmodel);
    1515void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads);
     16void solver_diagnostic_nonlinear_kgg(FemModel* femmodel,bool conserve_loads);
     17void solver_diagnostic_nonlinear_kff(FemModel* femmodel,bool conserve_loads);
    1618void solver_linear(FemModel* femmodel);
    1719void solver_adjoint_linear(FemModel* femmodel);

  • issm/trunk/src/m/classes/@model/model.m

    r5631 r5772  
    309309        md.qmu_relax=0;
    310310
     311        %solution parameters
     312        md.kff=0;
     313
    311314        %partitioner:
    312315        md.adjacency=[];

  • issm/trunk/src/m/classes/@model/setdefaultparameters.m

    r5631 r5772  
    256256%Ice solver: 'general' for Matlab's default solver (or 'lu' or 'sholesky')
    257257md.solver_type='general';
     258
     259%solution speed-up
     260md.kff=0;

  • issm/trunk/src/m/classes/public/marshall.m

    r5631 r5772  
    138138WriteData(fid,md.stokesreconditioning,'Scalar','stokesreconditioning');
    139139WriteData(fid,md.waitonlock,'Scalar','waitonlock');
     140WriteData(fid,md.kff,'Integer','kff');
    140141
    141142%Thermal parameters
Note: See TracChangeset for help on using the changeset viewer.