Index: ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 13876) +++ ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 13877) @@ -1,129 +0,0 @@ -/*!\file SystemMatricesx - * \brief: create system matrices (stiffness matrix, loads vector) - */ - -#include "./SystemMatricesx.h" -#include "../../shared/shared.h" -#include "../../include/include.h" -#include "../../io/io.h" -#include "../../toolkits/toolkits.h" -#include "../../EnumDefinitions/EnumDefinitions.h" - -void SystemMatricesx(Matrix** pKff, Matrix** pKfs, Vector** ppf, Vector** pdf, IssmDouble* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){ - - /*intermediary: */ - int i; - int fsize,ssize; - int connectivity, numberofdofspernode; - int analysis_type, configuration_type; - Element *element = NULL; - Load *load = NULL; - - /*output: */ - Matrix* Kff = NULL; - Matrix* Kfs = NULL; - Vector* pf = NULL; - Vector* df=NULL; - IssmDouble kmax = 0; - - /*Display message*/ - if(VerboseModule()) _pprintLine_(" Generating matrices"); - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum); - - /*Get size of matrices: */ - fsize=nodes->NumberOfDofs(configuration_type,FsetEnum); - ssize=nodes->NumberOfDofs(configuration_type,SsetEnum); - - numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum); - - /*Checks in debugging mode {{{*/ - if(penalty_kflag)_assert_(kflag); - if(penalty_pflag)_assert_(pflag); - /*}}}*/ - - /*Compute penalty free mstiffness matrix and load vector*/ - if(kflag){ - - Kff=new Matrix(fsize,fsize,connectivity,numberofdofspernode); - Kfs=new Matrix(fsize,ssize,connectivity,numberofdofspernode); - df=new Vector(fsize); - - /*Fill stiffness matrix from elements: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - element->CreateKMatrix(Kff,Kfs,df); - } - - /*Fill stiffness matrix from loads if loads have the current configuration_type: */ - for (i=0;iSize();i++){ - load=dynamic_cast(loads->GetObjectByOffset(i)); - if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs); - } - - /*Assemble matrix and doftypes and compress matrix to save memory: */ - Kff->Assemble(); - Kfs->Assemble(); - df->Assemble(); - } - - if(pflag){ - - pf=new Vector(fsize); - - /*Fill right hand side vector, from elements: */ - for (i=0;iSize();i++){ - element=dynamic_cast(elements->GetObjectByOffset(i)); - element->CreatePVector(pf); - } - - /*Fill right hand side from loads if loads have the current configuration_type: */ - for (i=0;iSize();i++){ - load=dynamic_cast(loads->GetObjectByOffset(i)); - if (load->InAnalysis(configuration_type)) load->CreatePVector(pf); - } - pf->Assemble(); - } - - /*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */ - kmax=Kff->Norm(NORM_INF); - - /*Now, deal with penalties*/ - if(penalty_kflag){ - - /*Fill stiffness matrix from loads: */ - for (i=0;iSize();i++){ - load=dynamic_cast(loads->GetObjectByOffset(i)); - if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax); - } - - /*Assemble matrix and compress matrix to save memory: */ - Kff->Assemble(); - Kfs->Assemble(); - } - - if(penalty_pflag){ - - /*Fill right hand side vector, from loads: */ - for (i=0;iSize();i++){ - load=dynamic_cast(loads->GetObjectByOffset(i)); - if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax); - } - - pf->Assemble(); - } - - /*Assign output pointers: */ - if(pKff) *pKff=Kff; - else xdelete(&Kff); - if(pKfs) *pKfs=Kfs; - else xdelete(&Kfs); - if(ppf) *ppf=pf; - else xdelete(&pf); - if(pdf) *pdf=df; - else xdelete(&df); - if(pkmax) *pkmax=kmax; -} Index: ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h =================================================================== --- ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h (revision 13876) +++ ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h (revision 13877) @@ -1,15 +0,0 @@ -/*!\file: SystemMatricesx.h - * \brief header file for degree of freedoms distribution routines. - */ - -#ifndef _SYSTEMMATRICESX_H -#define _SYSTEMMATRICESX_H - -#include "../../Container/Container.h" -#include "../../classes/objects/objects.h" - -/* local prototypes: */ -void SystemMatricesx(Matrix** pKff, Matrix** pKfs, Vector** ppf, Vector** pdf, IssmDouble* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters, - bool kflag=true,bool pflag=true,bool penalty_kflag=true,bool penalty_pflag=true); - -#endif /* _SYSTEMMATRICESX_H */ Index: ../trunk-jpl/src/c/modules/modules.h =================================================================== --- ../trunk-jpl/src/c/modules/modules.h (revision 13876) +++ ../trunk-jpl/src/c/modules/modules.h (revision 13877) @@ -90,7 +90,6 @@ #include "./Solverx/Solverx.h" #include "./SpcNodesx/SpcNodesx.h" #include "./SurfaceAreax/SurfaceAreax.h" -#include "./SystemMatricesx/SystemMatricesx.h" #include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h" #include "./TriaSearchx/TriaSearchx.h" #include "./TriMeshx/TriMeshx.h" Index: ../trunk-jpl/src/c/Makefile.am =================================================================== --- ../trunk-jpl/src/c/Makefile.am (revision 13876) +++ ../trunk-jpl/src/c/Makefile.am (revision 13877) @@ -298,8 +298,6 @@ ./modules/EnumToStringx/EnumToStringx.h\ ./modules/StringToEnumx/StringToEnumx.cpp\ ./modules/StringToEnumx/StringToEnumx.h\ - ./modules/SystemMatricesx/SystemMatricesx.cpp\ - ./modules/SystemMatricesx/SystemMatricesx.h\ ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp\ ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h\ ./modules/ConstraintsStatex/ConstraintsStatex.cpp\ Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 13876) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 13877) @@ -517,6 +517,110 @@ } } /*}}}*/ +void FemModel::SystemMatricesx(Matrix** pKff, Matrix** pKfs, Vector** ppf, Vector** pdf, IssmDouble* pkmax){/*{{{*/ + + /*intermediary: */ + int i; + int fsize,ssize; + int connectivity, numberofdofspernode; + int analysis_type, configuration_type; + Element *element = NULL; + Load *load = NULL; + + /*output: */ + Matrix *Kff = NULL; + Matrix *Kfs = NULL; + Vector *pf = NULL; + Vector *df = NULL; + IssmDouble kmax; + + /*Display message*/ + if(VerboseModule()) _pprintLine_(" Generating matrices"); + + /*retrive parameters: */ + this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + this->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum); + + /*Get size of matrices: */ + fsize=this->nodes->NumberOfDofs(configuration_type,FsetEnum); + ssize=this->nodes->NumberOfDofs(configuration_type,SsetEnum); + + numberofdofspernode=this->nodes->MaxNumDofs(configuration_type,GsetEnum); + + /*Compute penalty free mstiffness matrix and load vector*/ + Kff=new Matrix(fsize,fsize,connectivity,numberofdofspernode); + Kfs=new Matrix(fsize,ssize,connectivity,numberofdofspernode); + df =new Vector(fsize); + + /*Fill stiffness matrix from elements: */ + for (i=0;ielements->Size();i++){ + element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->CreateKMatrix(Kff,Kfs,df); + } + + /*Fill stiffness matrix from loads if loads have the current configuration_type: */ + for (i=0;iloads->Size();i++){ + load=dynamic_cast(this->loads->GetObjectByOffset(i)); + if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs); + } + + /*Assemble matrices*/ + Kff->Assemble(); + Kfs->Assemble(); + df->Assemble(); + + /*Create Load vector*/ + pf=new Vector(fsize); + + /*Fill right hand side vector, from elements: */ + for (i=0;ielements->Size();i++){ + element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->CreatePVector(pf); + } + + /*Fill right hand side from loads if loads have the current configuration_type: */ + for (i=0;iloads->Size();i++){ + load=dynamic_cast(this->loads->GetObjectByOffset(i)); + if (load->InAnalysis(configuration_type)) load->CreatePVector(pf); + } + pf->Assemble(); + + /*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */ + kmax=Kff->Norm(NORM_INF); + + /*Now, deal with penalties*/ + /*Fill stiffness matrix from loads: */ + for (i=0;iloads->Size();i++){ + load=dynamic_cast(this->loads->GetObjectByOffset(i)); + if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax); + } + + /*Assemble matrices*/ + Kff->Assemble(); + Kfs->Assemble(); + + /*Fill right hand side vector, from loads: */ + for (i=0;iloads->Size();i++){ + load=dynamic_cast(this->loads->GetObjectByOffset(i)); + if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax); + } + + /*Assemble vector*/ + pf->Assemble(); + + /*Assign output pointers: */ + if(pKff) *pKff=Kff; + else xdelete(&Kff); + if(pKfs) *pKfs=Kfs; + else xdelete(&Kfs); + if(ppf) *ppf=pf; + else xdelete(&pf); + if(pdf) *pdf=df; + else xdelete(&df); + if(pkmax) *pkmax=kmax; +} +/*}}}*/ void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/ int i; Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 13876) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 13877) @@ -88,6 +88,7 @@ void CostFunctionx( IssmDouble* pJ); void ThicknessAbsGradientx( IssmDouble* pJ, bool process_units,int weight_index); #endif + void SystemMatricesx(Matrix** pKff, Matrix** pKfs, Vector** ppf, Vector** pdf, IssmDouble* pkmax); void TimeAdaptx(IssmDouble* pdt); void UpdateConstraintsx(void); int UpdateVertexPositionsx(void); Index: ../trunk-jpl/src/c/solvers/solver_nonlinear.cpp =================================================================== --- ../trunk-jpl/src/c/solvers/solver_nonlinear.cpp (revision 13876) +++ ../trunk-jpl/src/c/solvers/solver_nonlinear.cpp (revision 13877) @@ -23,7 +23,7 @@ Vector* df = NULL; Vector* ys = NULL; - Loads* loads=NULL; + Loads* savedloads=NULL; bool converged; int constraints_converged; int num_unstable_constraints; @@ -32,7 +32,7 @@ /*parameters:*/ int min_mechanical_constraints; int max_nonlinear_iterations; - int configuration_type; + int configuration_type; /*Recover parameters: */ femmodel->parameters->FindParam(&min_mechanical_constraints,DiagnosticRiftPenaltyThresholdEnum); @@ -41,14 +41,15 @@ femmodel->UpdateConstraintsx(); /*Were loads requested as output? : */ - if(conserve_loads) loads=static_cast(femmodel->loads->Copy()); //protect loads from being modified by the solution - else loads=static_cast(femmodel->loads); //modify loads in this solution + if(conserve_loads){ + savedloads = static_cast(femmodel->loads->Copy()); + } count=1; converged=false; /*Start non-linear iteration using input velocity: */ - GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters); + GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters); //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) @@ -61,7 +62,7 @@ xdelete(&old_ug);old_ug=ug; xdelete(&old_uf);old_uf=uf; - SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters); @@ -71,7 +72,7 @@ InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); if(VerboseConvergence()) _pprintLine_(" number of unstable constraints: " << num_unstable_constraints); //rift convergence @@ -107,7 +108,10 @@ if(VerboseConvergence()) _pprintLine_("\n total number of iterations: " << count-1); /*clean-up*/ - if(conserve_loads) delete loads; + if(conserve_loads){ + delete femmodel->loads; + femmodel->loads=savedloads; + } xdelete(&uf); xdelete(&ug); xdelete(&old_ug); Index: ../trunk-jpl/src/c/solvers/solver_linear.cpp =================================================================== --- ../trunk-jpl/src/c/solvers/solver_linear.cpp (revision 13876) +++ ../trunk-jpl/src/c/solvers/solver_linear.cpp (revision 13877) @@ -23,7 +23,7 @@ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); femmodel->UpdateConstraintsx(); - SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); Index: ../trunk-jpl/src/c/solvers/solver_newton.cpp =================================================================== --- ../trunk-jpl/src/c/solvers/solver_newton.cpp (revision 13876) +++ ../trunk-jpl/src/c/solvers/solver_newton.cpp (revision 13877) @@ -55,7 +55,7 @@ xdelete(&old_uf);old_uf=uf; /*Solver forward model*/ - SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf,Kfs,ys);xdelete(&Kfs); Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);xdelete(&df); @@ -82,7 +82,7 @@ } /*Prepare next iteration using Newton's method*/ - SystemMatricesx(&Kff,&Kfs,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf,Kfs,ys); xdelete(&Kfs); Index: ../trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp =================================================================== --- ../trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp (revision 13876) +++ ../trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp (revision 13877) @@ -62,7 +62,7 @@ xdelete(&old_uf_horiz); old_uf_horiz=uf_horiz; /*solve: */ - SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf_horiz, Kfs_horiz, ys); xdelete(&Kfs_horiz); Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, df_horiz,femmodel->parameters); @@ -76,7 +76,7 @@ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); /*solve: */ - SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf_vert, Kfs_vert, ys); xdelete(&Kfs_vert); Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); xdelete(&Kff_vert); xdelete(&pf_vert); xdelete(&df_vert); Index: ../trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp =================================================================== --- ../trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp (revision 13876) +++ ../trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp (revision 13877) @@ -52,7 +52,7 @@ for(;;){ - SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); xdelete(&tf); Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters); Index: ../trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp =================================================================== --- ../trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp (revision 13876) +++ ../trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp (revision 13877) @@ -25,7 +25,7 @@ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); femmodel->UpdateConstraintsx(); - SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf, Kfs, ys,true); xdelete(&Kfs); //true means spc = 0 Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); xdelete(&Kff); xdelete(&pf); xdelete(&df);