Changeset 16126


Ignore:
Timestamp:
09/12/13 09:26:10 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moving SystemMatrices back to modules

Location:
issm/trunk-jpl/src/c
Files:
9 added
12 edited

Legend:

Unmodified
Added
Removed

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

    r16113 r16126  
    279279                                        ./modules/SurfaceAreax/SurfaceAreax.h\
    280280                                        ./modules/SurfaceAreax/SurfaceAreax.cpp\
     281                                        ./modules/AllocateSystemMatricesx/AllocateSystemMatricesx.h\
     282                                        ./modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp\
     283                                        ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h\
     284                                        ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp\
     285                                        ./modules/SystemMatricesx/SystemMatricesx.h\
     286                                        ./modules/SystemMatricesx/SystemMatricesx.cpp\
    281287                                        ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h\
    282288                                        ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp\

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r16120 r16126  
    366366/*}}}*/
    367367/*Modules:*/
    368 void FemModel::AllocateSystemMatrices(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf){ /*{{{*/
    369 
    370         /*Intermediary*/
    371         int  fsize,ssize,flocalsize,slocalsize;
    372         int  connectivity, numberofdofspernode;
    373         int  configuration_type;
    374         int  m,n,M,N;
    375         int *d_nnz = NULL;
    376         int *o_nnz = NULL;
    377 
    378         /*output*/
    379         Matrix<IssmDouble> *Kff  = NULL;
    380         Matrix<IssmDouble> *Kfs  = NULL;
    381         Vector<IssmDouble> *pf   = NULL;
    382         Vector<IssmDouble> *df   = NULL;
    383 
    384         bool oldalloc=false;
    385 
    386         /*retrieve parameters: */
    387         this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    388         this->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
    389 
    390         /*retrieve node info*/
    391         fsize      = this->nodes->NumberOfDofs(configuration_type,FsetEnum);
    392         ssize      = this->nodes->NumberOfDofs(configuration_type,SsetEnum);
    393         flocalsize = this->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
    394         slocalsize = this->nodes->NumberOfDofsLocal(configuration_type,SsetEnum);
    395 
    396         numberofdofspernode=this->nodes->MaxNumDofs(configuration_type,GsetEnum);
    397 
    398         if(oldalloc){
    399                 if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
    400                 if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
    401                 if(pdf)  df =new Vector<IssmDouble>(fsize);
    402                 if(ppf)  pf =new Vector<IssmDouble>(fsize);
    403         }
    404         else{
    405                 if(pKff){
    406                         m=flocalsize; n=flocalsize; /*local  sizes*/
    407                         M=fsize;      N=fsize;      /*global sizes*/
    408                         this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,FsetEnum);
    409                         Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
    410                         xDelete<int>(d_nnz);
    411                         xDelete<int>(o_nnz);
    412                 }
    413                 if(pKfs){
    414                         m=flocalsize; n=slocalsize; /*local  sizes*/
    415                         M=fsize;      N=ssize;      /*global sizes*/
    416                         this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,SsetEnum);
    417                         Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
    418                         xDelete<int>(d_nnz);
    419                         xDelete<int>(o_nnz);
    420                 }
    421                 if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
    422                 if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
    423         }
    424 
    425         /*Allocate output pointers*/
    426         if(pKff) *pKff = Kff;
    427         if(pKfs) *pKfs = Kfs;
    428         if(pdf)  *pdf  = df;
    429         if(ppf)  *ppf  = pf;
    430 }
    431 /*}}}*/
    432 void FemModel::MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum){/*{{{*/
    433 
    434         /*Intermediary*/
    435         int      i,j,k,index,offset,count;
    436         int      configuration_type;
    437         int      d_nz,o_nz;
    438         Element *element            = NULL;
    439         Load    *load               = NULL;
    440         int     *head_e             = NULL;
    441         int     *next_e             = NULL;
    442         int     *count2offset_e     = NULL;
    443         int     *head_l             = NULL;
    444         int     *next_l             = NULL;
    445         int     *count2offset_l     = NULL;
    446         int     *lidlist            = NULL;
    447 
    448         /*output*/
    449         int *d_nnz = NULL;
    450         int *o_nnz = NULL;
    451 
    452         /*retrive parameters: */
    453         this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    454 
    455         /*Get vector size and number of nodes*/
    456         int numnodes            = nodes->NumberOfNodes(configuration_type);
    457         int localnumnodes       = nodes->Size();
    458         int numberofdofspernode = nodes->MaxNumDofs(configuration_type,GsetEnum);
    459         int M                   = nodes->NumberOfDofs(configuration_type,set1enum);
    460         int N                   = nodes->NumberOfDofs(configuration_type,set2enum);
    461         int m                   = nodes->NumberOfDofsLocal(configuration_type,set1enum);
    462         int n                   = nodes->NumberOfDofsLocal(configuration_type,set2enum);
    463         int numnodesperelement  = elements->MaxNumNodes();
    464         int numnodesperload     = loads->MaxNumNodes(configuration_type);
    465 
    466         /*First, we are building chaining vectors so that we know what nodes are
    467          * connected to what elements. These vectors are such that:
    468          *   for(int i=head[id];i!=-1;i=next[i])
    469          * will loop over all the elements that are connected to the node number
    470          * id*/
    471         head_e         = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
    472         next_e         = xNew<int>(elements->Size()*numnodesperelement);
    473         count2offset_e = xNew<int>(elements->Size()*numnodesperelement);
    474 
    475         k=0;
    476         for(i=0;i<elements->Size();i++){
    477                 element = dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    478                 lidlist = xNew<int>(element->GetNumberOfNodes());
    479                 element->GetNodesLidList(lidlist);
    480 
    481                 for(j=0;j<element->GetNumberOfNodes();j++){
    482                         index = lidlist[j];
    483                         _assert_(index>=0 && index<numnodes);
    484 
    485                         count2offset_e[k]=i;
    486                         next_e[k]=head_e[index];
    487                         head_e[index]=k++;
    488                 }
    489                 for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;
    490 
    491                 xDelete<int>(lidlist);
    492         }
    493 
    494         /*Chain for loads*/
    495         head_l         = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
    496         next_l         = xNew<int>(loads->Size(configuration_type)*numnodesperload);
    497         count2offset_l = xNew<int>(loads->Size(configuration_type)*numnodesperload);
    498         k=0;
    499         for(i=0;i<loads->Size();i++){
    500                 load = dynamic_cast<Load*>(loads->GetObjectByOffset(i));
    501                 if(!load->InAnalysis(configuration_type)) continue;
    502                 lidlist = xNew<int>(load->GetNumberOfNodes());
    503                 load->GetNodesLidList(lidlist);
    504 
    505                 for(j=0;j<load->GetNumberOfNodes();j++){
    506                         index = lidlist[j];
    507                         _assert_(index>=0 && index<numnodes);
    508 
    509                         count2offset_l[k]=i;
    510                         next_l[k]=head_l[index];
    511                         head_l[index]=k++;
    512                 }
    513                 for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;
    514 
    515                 xDelete<int>(lidlist);
    516         }
    517 
    518         /*OK now count number of dofs and flag each nodes for each node i*/
    519         bool *flags                  = xNew<bool>(localnumnodes);
    520         int  *flagsindices           = xNew<int>(localnumnodes);
    521         int  *d_connectivity         = xNewZeroInit<int>(numnodes);
    522         int  *o_connectivity         = xNewZeroInit<int>(numnodes);
    523         int  *connectivity_clone     = xNewZeroInit<int>(numnodes);
    524         int  *all_connectivity_clone = xNewZeroInit<int>(numnodes);
    525 
    526         /*Resetting flags to false at eahc iteration takes a lot of time, so we keep track of the flags
    527          * to reset in flagsindices, initialized with -1*/
    528         for(i = 0;i<localnumnodes;i++) flags[i]        = false;
    529         for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
    530 
    531         /*Create connectivity vector*/
    532         for(i=0;i<nodes->Size();i++){
    533                 Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
    534                 if(node->InAnalysis(configuration_type)){
    535 
    536                         /*Reinitialize flags to false*/
    537                         j=0;
    538                         while(true){
    539                                 if(flagsindices[j]>=0){
    540                                         flags[flagsindices[j]] = false;
    541                                         flagsindices[j]        = -1;
    542                                         j++;
    543                                 }
    544                                 else{
    545                                         break;
    546                                 }
    547                         }
    548 
    549                         //for(j=0;j<localnumnodes;j++) flags[j]=false;
    550 
    551                         /*Loop over elements that hold node number i*/
    552                         //if(head_e[node->Lid()]==-1 && head_l[node->Lid()]==-1){
    553                         //      printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Lid()+1);
    554                         //}
    555                         for(j=head_e[node->Lid()];j!=-1;j=next_e[j]){
    556                                 offset=count2offset_e[j];
    557                                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(offset));
    558                                 element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
    559                                 if(node->IsClone()){
    560                                         connectivity_clone[node->Sid()]+=d_nz+o_nz;
    561                                 }
    562                                 else{
    563                                         d_connectivity[node->Sid()]+=d_nz;
    564                                         o_connectivity[node->Sid()]+=o_nz;
    565                                 }
    566                         }
    567                         for(j=head_l[node->Lid()];j!=-1;j=next_l[j]){
    568                                 offset=count2offset_l[j];
    569                                 load=dynamic_cast<Load*>(loads->GetObjectByOffset(offset));
    570                                 load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
    571                                 if(node->IsClone()){
    572                                         connectivity_clone[node->Sid()]+=d_nz+o_nz;
    573                                 }
    574                                 else{
    575                                         d_connectivity[node->Sid()]+=d_nz;
    576                                         o_connectivity[node->Sid()]+=o_nz;
    577                                 }
    578                         }
    579                 }
    580         }
    581         xDelete<bool>(flags);
    582         xDelete<int>(flagsindices);
    583         xDelete<int>(count2offset_e);
    584         xDelete<int>(head_e);
    585         xDelete<int>(next_e);
    586         xDelete<int>(count2offset_l);
    587         xDelete<int>(head_l);
    588         xDelete<int>(next_l);
    589 
    590         /*sum over all cpus*/
    591         ISSM_MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
    592         xDelete<int>(connectivity_clone);
    593 
    594         if(set1enum==FsetEnum){
    595                 count=0;
    596                 d_nnz=xNew<int>(m);
    597                 o_nnz=xNew<int>(m);
    598                 for(i=0;i<nodes->Size();i++){
    599                         Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
    600                         if(node->InAnalysis(configuration_type) && !node->IsClone()){
    601                                 for(j=0;j<node->indexing.fsize;j++){
    602                                         _assert_(count<m);
    603                                         d_nnz[count]=numberofdofspernode*(d_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
    604                                         o_nnz[count]=numberofdofspernode*(o_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
    605                                         if(d_nnz[count]>n)   d_nnz[count]=n;
    606                                         if(o_nnz[count]>N-n) o_nnz[count]=N-n;
    607                                         count++;
    608                                 }
    609                         }
    610                 }
    611                 _assert_(m==count);
    612         }
    613         else{
    614                 _error_("STOP not implemented");
    615         }
    616         xDelete<int>(d_connectivity);
    617         xDelete<int>(o_connectivity);
    618         xDelete<int>(all_connectivity_clone);
    619 
    620         /*Allocate ouptput pointer*/
    621         *pd_nnz=d_nnz;
    622         *po_nnz=o_nnz;
    623 
    624 }/*}}}*/
    625 void FemModel::CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,IssmDouble kmax){/*{{{*/
    626 
    627         int      i,connectivity;
    628         int      numberofdofspernode;
    629         int      fsize,configuration_type;
    630         Element *element = NULL;
    631         Load    *load    = NULL;
    632         Matrix<IssmDouble>* Jff = NULL;
    633 
    634         /*Checks*/
    635         _assert_(nodes && elements);
    636 
    637         /*Recover some parameters*/
    638         parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    639         parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
    640         fsize=nodes->NumberOfDofs(configuration_type,FsetEnum);
    641         numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum);
    642 
    643         /*Initialize Jacobian Matrix*/
    644         this->AllocateSystemMatrices(&Jff,NULL,NULL,NULL);
    645 
    646         /*Create and assemble matrix*/
    647         for(i=0;i<elements->Size();i++){
    648                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    649                 element->CreateJacobianMatrix(Jff);
    650         }
    651         for (i=0;i<loads->Size();i++){
    652                 load=(Load*)loads->GetObjectByOffset(i);
    653                 if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff);
    654                 if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax);
    655         }
    656         Jff->Assemble();
    657 
    658         /*Assign output pointer*/
    659         *pJff=Jff;
    660 
    661 }/*}}}*/
    662368int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
    663369
    … …  
    845551}
    846552/*}}}*/
    847 void FemModel::SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax){/*{{{*/
    848 
    849         /*intermediary: */
    850         int      i,M,N;
    851         int      configuration_type;
    852         Element *element = NULL;
    853         Load    *load    = NULL;
    854 
    855         /*output: */
    856         Matrix<IssmDouble> *Kff  = NULL;
    857         Matrix<IssmDouble> *Kfs  = NULL;
    858         Vector<IssmDouble> *pf   = NULL;
    859         Vector<IssmDouble> *df   = NULL;
    860         IssmDouble          kmax = 0;
    861 
    862         /*Display message*/
    863         if(VerboseModule()) _printf0_("   Generating matrices");
    864 
    865         /*retrive parameters: */
    866         this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    867 
    868         /*First, we might need to do a dry run to get kmax if penalties are employed*/
    869         if(loads->IsPenalty(configuration_type)){
    870 
    871                 /*Allocate Kff_temp*/
    872                 Matrix<IssmDouble> *Kff_temp = NULL;
    873                 this->AllocateSystemMatrices(&Kff_temp,NULL,NULL,NULL);
    874 
    875                 /*Get complete stiffness matrix without penalties*/
    876                 for (i=0;i<this->elements->Size();i++){
    877                         element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    878                         element->CreateKMatrix(Kff_temp,NULL);
    879                 }
    880 
    881                 for (i=0;i<this->loads->Size();i++){
    882                         load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    883                         if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL);
    884                 }
    885                 Kff_temp->Assemble();
    886 
    887                 /*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */
    888                 kmax=Kff_temp->Norm(NORM_INF);
    889                 delete Kff_temp;
    890         }
    891 
    892         /*Allocate stiffness matrices and load vector*/
    893         this->AllocateSystemMatrices(&Kff,&Kfs,&df,&pf);
    894 
    895         /*Display size*/
    896         if(VerboseModule()){
    897                 Kff->GetSize(&M,&N);
    898                 _printf0_(" (Kff stiffness matrix size: "<<M<<" x "<<N<<")\n");
    899         }
    900 
    901         /*Fill stiffness matrix from elements and loads */
    902         for (i=0;i<this->elements->Size();i++){
    903                 element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    904                 element->CreateKMatrix(Kff,Kfs);
    905         }
    906 
    907         for (i=0;i<this->loads->Size();i++){
    908                 load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    909                 if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
    910         }
    911 
    912         /*Fill right hand side vector, from elements and loads */
    913         for (i=0;i<this->elements->Size();i++){
    914                 element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    915                 element->CreatePVector(pf);
    916         }
    917         for (i=0;i<this->loads->Size();i++){
    918                 load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    919                 if(load->InAnalysis(configuration_type)) load->CreatePVector(pf);
    920         }
    921 
    922         /*Now deal with penalties (only in loads)*/
    923         if(loads->IsPenalty(configuration_type)){
    924                 for (i=0;i<this->loads->Size();i++){
    925                         load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    926                         if(load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
    927                 }
    928                 for (i=0;i<this->loads->Size();i++){
    929                         load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
    930                         if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
    931                 }
    932         }
    933 
    934         /*Create dof vector for stiffness matrix preconditioning*/
    935         for (i=0;i<this->elements->Size();i++){
    936                 element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    937                 element->CreateDVector(df);
    938         }
    939 
    940         /*Assemble matrices and vector*/
    941         Kff->Assemble();
    942         Kfs->Assemble();
    943         pf->Assemble();
    944         df->Assemble();
    945         //Kff->AllocationInfo();
    946         //Kfs->AllocationInfo();
    947 
    948         /*Assign output pointers: */
    949         if(pKff) *pKff=Kff;
    950         else      delete Kff;
    951         if(pKfs) *pKfs=Kfs;
    952         else      delete Kfs;
    953         if(ppf)  *ppf=pf;
    954         else      delete pf;
    955         if(pdf)  *pdf=df;
    956         else      delete df;
    957         if(pkmax) *pkmax=kmax;
    958 }
    959 /*}}}*/
    960553void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
    961554

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r15853 r16126  
    5151
    5252                /*Methods:*/
    53                 void AllocateSystemMatrices(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf);
    54                 void CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,IssmDouble kmax);
    5553                void Echo();
    5654                void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
    57                 void MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum);
    5855                void Solve(void);
    5956                void OutputResults(void);
    … …  
    9592                void Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y);
    9693                #endif
    97                 void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax);
    9894                void TimeAdaptx(IssmDouble* pdt);
    9995                void UpdateConstraintsx(void);

  • issm/trunk-jpl/src/c/modules/Exp2Kmlx/Exp2Kmlx.h

    r15000 r16126  
    1010
    1111/* local prototypes: */
    12 int Exp2Kmlx(char* filexp,char* filkml,
    13                          int sgn,
    14                          bool holes);
    15 int Exp2Kmlx(char* filexp,char* filkml,
    16                          int sgn,double cm,double sp,
    17                          bool holes);
     12int Exp2Kmlx(char* filexp,char* filkml,int sgn,bool holes);
     13int Exp2Kmlx(char* filexp,char* filkml,int sgn,double cm,double sp,bool holes);
    1814
    1915#endif  /* _EXP2KMLX_H */

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

    r16026 r16126  
    77
    88/*Modules: */
     9#include "./AllocateSystemMatricesx/AllocateSystemMatricesx.h"
    910#include "./AverageFilterx/AverageFilterx.h"
    1011#include "./AverageOntoPartitionx/AverageOntoPartitionx.h"
    … …  
    2223#include "./ControlInputScaleGradientx/ControlInputScaleGradientx.h"
    2324#include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h"
     25#include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
    2426#include "./Delta18oParameterizationx/Delta18oParameterizationx.h"
    2527#include "./DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
    … …  
    8890#include "./SmbGradientsx/SmbGradientsx.h"
    8991#include "./Solverx/Solverx.h"
     92#include "./SystemMatricesx/SystemMatricesx.h"
    9093#include "./SpcNodesx/SpcNodesx.h"
    9194#include "./SurfaceAreax/SurfaceAreax.h"

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_adjoint_linear.cpp

    r15849 r16126  
    2626        femmodel->UpdateConstraintsx();
    2727
    28         femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
     28        SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel);
    2929        CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    3030        Reduceloadx(pf, Kfs, ys,true); delete Kfs; //true means spc = 0

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

    r15849 r16126  
    7979                for(;;){
    8080                        femmodel->HydrologyTransferx();
    81                         femmodel->SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax);
     81                        SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel);
    8282                        CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
    8383                        Reduceloadx(pf,Kfs,ys); delete Kfs;
    … …  
    123123                        for(;;){
    124124                                femmodel->HydrologyTransferx();
    125                                 femmodel->SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL);
     125                                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    126126                                CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
    127127                                Reduceloadx(pf,Kfs,ys); delete Kfs;

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp

    r15849 r16126  
    2424        femmodel->UpdateConstraintsx();
    2525
    26         femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
     26        SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    2727        CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    2828        Reduceloadx(pf, Kfs, ys); delete Kfs;

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp

    r15849 r16126  
    5656                /*Solver forward model*/
    5757                if(count==1 || newton==2){
    58                         femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
     58                        SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    5959                        CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    6060                        Reduceloadx(pf,Kfs,ys);delete Kfs;
    … …  
    6868
    6969                /*Prepare next iteration using Newton's method*/
    70                 femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);delete df;
     70                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);delete df;
    7171                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    7272                Reduceloadx(pf,Kfs,ys);delete Kfs;
    … …  
    7676                pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);
    7777
    78                 femmodel->CreateJacobianMatrixx(&Jff,kmax);
     78                CreateJacobianMatrixx(&Jff,femmodel,kmax);
    7979                Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); delete Jff; delete pJf;
    8080                uf->AXPY(duf, 1.0); delete duf;

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp

    r15849 r16126  
    6060                delete ug;
    6161
    62                 femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
     62                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    6363                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    6464                Reduceloadx(pf, Kfs, ys); delete Kfs;

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp

    r15849 r16126  
    6161
    6262                /*solve: */
    63                 femmodel->SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL);
     63                SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL,femmodel);
    6464                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    6565                Reduceloadx(pf_horiz, Kfs_horiz, ys); delete Kfs_horiz;
    … …  
    7575
    7676                /*solve: */
    77                 femmodel->SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert,  &df_vert,NULL);
     77                SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL,femmodel);
    7878                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    7979                Reduceloadx(pf_vert, Kfs_vert, ys); delete Kfs_vert;

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp

    r15849 r16126  
    4848
    4949                delete tf_old; tf_old=tf;
    50                 femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset);
     50                SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
    5151                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
    5252                Reduceloadx(pf, Kfs, ys); delete Kfs;
Note: See TracChangeset for help on using the changeset viewer.