Index: /issm/trunk-jpl/src/c/Container/Loads.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Loads.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/Container/Loads.cpp (revision 13925) @@ -51,4 +51,31 @@ } +} +/*}}}*/ +/*FUNCTION Loads::IsPenalty{{{*/ +bool Loads::IsPenalty(int analysis_type){ + + int ispenalty=0; + int allispenalty=0; + + /*Now go through all loads, and get how many nodes they own, unless they are clone nodes: */ + for(int i=0;iSize();i++){ + + Load* load=dynamic_cast(this->GetObjectByOffset(i)); + if (load->InAnalysis(analysis_type)){ + if(load->IsPenalty()) ispenalty++; + } + } + + /*Grab sum of all cpus: */ +#ifdef _HAVE_MPI_ + MPI_Allreduce((void*)&ispenalty,(void*)&allispenalty,1,MPI_INT,MPI_SUM,IssmComm::GetComm()); + ispenalty=allispenalty; +#endif + + if(ispenalty) + return true; + else + return false; } /*}}}*/ Index: /issm/trunk-jpl/src/c/Container/Loads.h =================================================================== --- /issm/trunk-jpl/src/c/Container/Loads.h (revision 13924) +++ /issm/trunk-jpl/src/c/Container/Loads.h (revision 13925) @@ -26,4 +26,5 @@ /*numerics*/ void Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); + bool IsPenalty(int analysis); int MaxNumNodes(int analysis); int NumberOfLoads(void); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 13925) @@ -353,5 +353,4 @@ Vector *df = NULL; - //bool oldalloc=false; bool oldalloc=true; @@ -369,36 +368,35 @@ if(oldalloc){ - Kff=new Matrix(fsize,fsize,connectivity,numberofdofspernode); - Kfs=new Matrix(fsize,ssize,connectivity,numberofdofspernode); - df =new Vector(fsize); - pf =new Vector(fsize); + if(pKff) Kff=new Matrix(fsize,fsize,connectivity,numberofdofspernode); + if(pKfs) Kfs=new Matrix(fsize,ssize,connectivity,numberofdofspernode); + if(pdf) df =new Vector(fsize); + if(ppf) pf =new Vector(fsize); } else{ - /*Kff*/ - m=flocalsize; n=flocalsize; /*local sizes*/ - M=fsize; N=fsize; /*global sizes*/ - this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,FsetEnum); - Kff=new Matrix(m,n,M,N,d_nnz,o_nnz); - xDelete(d_nnz); - xDelete(o_nnz); - - /*Kfs*/ - m=flocalsize; n=slocalsize; /*local sizes*/ - M=fsize; N=ssize; /*global sizes*/ - this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,SsetEnum); - Kfs=new Matrix(m,n,M,N,d_nnz,o_nnz); - xDelete(d_nnz); - xDelete(o_nnz); - - /*Vectors*/ - df =new Vector(flocalsize,fsize); - pf =new Vector(flocalsize,fsize); + if(pKff){ + m=flocalsize; n=flocalsize; /*local sizes*/ + M=fsize; N=fsize; /*global sizes*/ + this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,FsetEnum); + Kff=new Matrix(m,n,M,N,d_nnz,o_nnz); + xDelete(d_nnz); + xDelete(o_nnz); + } + if(pKfs){ + m=flocalsize; n=slocalsize; /*local sizes*/ + M=fsize; N=ssize; /*global sizes*/ + this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,SsetEnum); + Kfs=new Matrix(m,n,M,N,d_nnz,o_nnz); + xDelete(d_nnz); + xDelete(o_nnz); + } + if(pdf) df =new Vector(flocalsize,fsize); + if(ppf) pf =new Vector(flocalsize,fsize); } /*Allocate output pointers*/ - *pKff = Kff; - *pKfs = Kfs; - *pdf = df; - *ppf = pf; + if(pKff) *pKff = Kff; + if(pKfs) *pKfs = Kfs; + if(pdf) *pdf = df; + if(ppf) *ppf = pf; } /*}}}*/ @@ -770,5 +768,5 @@ Vector *pf = NULL; Vector *df = NULL; - IssmDouble kmax; + IssmDouble kmax = 0; /*Display message*/ @@ -778,5 +776,28 @@ this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - /*Compute penalty free mstiffness matrix and load vector*/ + /*First, we might need to do a dry run to get kmax if penalties are employed*/ + if(loads->IsPenalty(configuration_type)){ + + /*Allocate Kff_temp*/ + Matrix *Kff_temp = NULL; + this->AllocateSystemMatrices(&Kff_temp,NULL,NULL,NULL); + + /*Get complete stiffness matrix without penalties*/ + for (i=0;ielements->Size();i++){ + element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->CreateKMatrix(Kff_temp,NULL,NULL); + } + for (i=0;iloads->Size();i++){ + load=dynamic_cast(this->loads->GetObjectByOffset(i)); + if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL); + } + Kff_temp->Assemble(); + + /*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */ + kmax=Kff_temp->Norm(NORM_INF); + delete Kff_temp; + } + + /*Allocate stiffness matrices and load vector*/ this->AllocateSystemMatrices(&Kff,&Kfs,&df,&pf); @@ -788,5 +809,5 @@ for (i=0;iloads->Size();i++){ load=dynamic_cast(this->loads->GetObjectByOffset(i)); - if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs); + if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs); } @@ -798,23 +819,17 @@ for (i=0;iloads->Size();i++){ load=dynamic_cast(this->loads->GetObjectByOffset(i)); - if (load->InAnalysis(configuration_type)) load->CreatePVector(pf); - } - - /*Assemble matrices (required to get kmax)*/ - Kff->Assemble(); - Kfs->Assemble(); - df->Assemble(); - - /*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */ - kmax=Kff->Norm(NORM_INF); - - /*Now that we have kmax, deal with penalties (only in loads)*/ - for (i=0;iloads->Size();i++){ - load=dynamic_cast(this->loads->GetObjectByOffset(i)); - if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax); - } - for (i=0;iloads->Size();i++){ - load=dynamic_cast(this->loads->GetObjectByOffset(i)); - if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax); + if(load->InAnalysis(configuration_type)) load->CreatePVector(pf); + } + + /*Now deal with penalties (only in loads)*/ + if(loads->IsPenalty(configuration_type)){ + for (i=0;iloads->Size();i++){ + load=dynamic_cast(this->loads->GetObjectByOffset(i)); + if(load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax); + } + for (i=0;iloads->Size();i++){ + load=dynamic_cast(this->loads->GetObjectByOffset(i)); + if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax); + } } @@ -823,4 +838,5 @@ Kfs->Assemble(); pf->Assemble(); + df->Assemble(); //Kff->AllocationInfo(); //Kfs->AllocationInfo(); Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp (revision 13925) @@ -316,4 +316,9 @@ } +} +/*}}}*/ +/*FUNCTION Icefront::IsPenalty{{{*/ +bool Icefront::IsPenalty(void){ + return false; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.h (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.h (revision 13925) @@ -73,4 +73,5 @@ int GetNumberOfNodes(void); void GetNodesSidList(int* sidlist); + bool IsPenalty(void); void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, IssmDouble kmax); void PenaltyCreatePVector(Vector* pf, IssmDouble kmax); Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Load.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Load.h (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Load.h (revision 13925) @@ -25,4 +25,5 @@ virtual ~Load(){}; virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; + virtual bool IsPenalty(void)=0; virtual int GetNumberOfNodes(void)=0; virtual void GetNodesSidList(int* sidlist)=0; Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.cpp (revision 13925) @@ -340,4 +340,9 @@ } +} +/*}}}*/ +/*FUNCTION Numericalflux::IsPenalty{{{*/ +bool Numericalflux::IsPenalty(void){ + return false; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.h (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.h (revision 13925) @@ -67,4 +67,5 @@ int GetNumberOfNodes(void); void CreateJacobianMatrix(Matrix* Jff){_error_("Not implemented yet");}; + bool IsPenalty(void); void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax){_error_("Not implemented yet");}; void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, IssmDouble kmax); Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp (revision 13925) @@ -293,4 +293,9 @@ if (in_analysis_type==this->analysis_type)return true; else return false; +} +/*}}}*/ +/*FUNCTION Pengrid::IsPenalty{{{*/ +bool Pengrid::IsPenalty(void){ + return true; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.h (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.h (revision 13925) @@ -74,4 +74,5 @@ void GetNodesSidList(int* sidlist); int GetNumberOfNodes(void); + bool IsPenalty(void); void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax){_error_("Not implemented yet");}; void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, IssmDouble kmax); Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.cpp (revision 13925) @@ -161,4 +161,9 @@ return NUMVERTICES; +} +/*}}}*/ +/*FUNCTION Penpair::IsPenalty{{{*/ +bool Penpair::IsPenalty(void){ + return true; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.h (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.h (revision 13925) @@ -59,4 +59,5 @@ void GetNodesSidList(int* sidlist); int GetNumberOfNodes(void); + bool IsPenalty(void); void PenaltyCreateKMatrix(Matrix* Kff,Matrix* Kfs,IssmDouble kmax); void PenaltyCreatePVector(Vector* pf, IssmDouble kmax); Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.cpp (revision 13925) @@ -293,4 +293,9 @@ this->parameters=parametersin; +} +/*}}}*/ +/*FUNCTION Riftfront::IsPenalty{{{*/ +bool Riftfront::IsPenalty(void){ + return true; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.h (revision 13924) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.h (revision 13925) @@ -80,4 +80,5 @@ void GetNodesSidList(int* sidlist); int GetNumberOfNodes(void); + bool IsPenalty(void); void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax){_error_("Not implemented yet");}; void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, IssmDouble kmax); Index: /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp (revision 13925) @@ -12,12 +12,16 @@ /*variables: */ - int i; - int configuration_type; - int fsize; - IssmDouble* ug_serial=NULL; + int configuration_type; + int fsize; + IssmDouble *ug_serial = NULL; + bool oldalloc = true; /*first figure out fsize: */ parameters->FindParam(&configuration_type,ConfigurationTypeEnum); fsize=nodes->NumberOfDofs(configuration_type,FsetEnum); + + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + int flocalsize = nodes->NumberOfDofsLocal(analysis_type,FsetEnum); if(fsize==0){ @@ -26,5 +30,8 @@ else{ /*allocate: */ - uf=new Vector(fsize); + if(oldalloc) + uf=new Vector(fsize); + else + uf=new Vector(flocalsize,fsize); if(nodes->NumberOfNodes(configuration_type)){ @@ -34,5 +41,5 @@ /*Go through all nodes, and ask them to retrieve values from ug, and plug them into uf: */ - for(i=0;iSize();i++){ + for(int i=0;iSize();i++){ Node* node=(Node*)nodes->GetObjectByOffset(i); Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp (revision 13924) +++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp (revision 13925) @@ -40,5 +40,5 @@ } /*}}}*/ -/*FUNCTION PetscMat::PetscMat(int M,int N, IssmDouble sparsity){{{*/ +/*FUNCTION PetscMat::PetscMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/ PetscMat::PetscMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){