Changeset 13877

Timestamp:

11/05/12 13:56:10 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moved SystemMatrices to FemModel.cpp (for solver_nonlinear, needed to change femmodel->load if conserveloads=true)

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (1 diff)
• classes/FemModel.cpp (modified) (1 diff)
• classes/FemModel.h (modified) (1 diff)
• modules/SystemMatricesx (deleted)
• modules/modules.h (modified) (1 diff)
• solvers/solver_adjoint_linear.cpp (modified) (1 diff)
• solvers/solver_linear.cpp (modified) (1 diff)
• solvers/solver_newton.cpp (modified) (2 diffs)
• solvers/solver_nonlinear.cpp (modified) (7 diffs)
• solvers/solver_stokescoupling_nonlinear.cpp (modified) (2 diffs)
• solvers/solver_thermal_nonlinear.cpp (modified) (1 diff)

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

@@ -299 +299 @@
                                         ./modules/StringToEnumx/StringToEnumx.cpp\
                                         ./modules/StringToEnumx/StringToEnumx.h\
-                                        ./modules/SystemMatricesx/SystemMatricesx.cpp\
-                                        ./modules/SystemMatricesx/SystemMatricesx.h\
                                         ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp\
                                         ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h\

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

@@ -518 +518 @@
 }
 /*}}}*/
+void FemModel::SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** 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<IssmDouble> *Kff  = NULL;
+        Matrix<IssmDouble> *Kfs  = NULL;
+        Vector<IssmDouble> *pf   = NULL;
+        Vector<IssmDouble> *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<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
+        Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
+        df =new Vector<IssmDouble>(fsize);
+
+        /*Fill stiffness matrix from elements: */
+        for (i=0;i<this->elements->Size();i++){
+                element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                element->CreateKMatrix(Kff,Kfs,df);
+        }
+
+        /*Fill stiffness matrix from loads if loads have the current configuration_type: */
+        for (i=0;i<this->loads->Size();i++){
+                load=dynamic_cast<Load*>(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<IssmDouble>(fsize);
+
+        /*Fill right hand side vector, from elements: */
+        for (i=0;i<this->elements->Size();i++){
+                element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                element->CreatePVector(pf);
+        }
+
+        /*Fill right hand side from loads if loads have the current configuration_type: */
+        for (i=0;i<this->loads->Size();i++){
+                load=dynamic_cast<Load*>(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;i<this->loads->Size();i++){
+                load=dynamic_cast<Load*>(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;i<this->loads->Size();i++){
+                load=dynamic_cast<Load*>(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){/*{{{*/
 

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

@@ -89 +89 @@
                 void ThicknessAbsGradientx( IssmDouble* pJ, bool process_units,int weight_index);
                 #endif
+                void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax);
                 void TimeAdaptx(IssmDouble* pdt);
                 void UpdateConstraintsx(void);

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

@@ -91 +91 @@
 #include "./SpcNodesx/SpcNodesx.h"
 #include "./SurfaceAreax/SurfaceAreax.h"
-#include "./SystemMatricesx/SystemMatricesx.h"
 #include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
 #include "./TriaSearchx/TriaSearchx.h"

issm/trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp

@@ -26 +26 @@
         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

issm/trunk-jpl/src/c/solvers/solver_linear.cpp

@@ -24 +24 @@
         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);

issm/trunk-jpl/src/c/solvers/solver_newton.cpp

@@ -56 +56 @@
 
                 /*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);
@@ -83 +83 @@
 
                 /*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);

issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp

@@ -24 +24 @@
         Vector<IssmDouble>* ys  = NULL;
 
-        Loads* loads=NULL;
+        Loads* savedloads=NULL;
         bool converged;
         int constraints_converged;
@@ -33 +33 @@
         int min_mechanical_constraints;
         int max_nonlinear_iterations;
-        int  configuration_type;
+        int configuration_type;
 
         /*Recover parameters: */
@@ -42 +42 @@
 
         /*Were loads requested as output? : */
-        if(conserve_loads) loads=static_cast<Loads*>(femmodel->loads->Copy()); //protect loads from being modified by the solution
-        else               loads=static_cast<Loads*>(femmodel->loads);         //modify loads  in this solution
+        if(conserve_loads){
+                savedloads = static_cast<Loads*>(femmodel->loads->Copy());
+        }
 
         count=1;
@@ -49 +50 @@
 
         /*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);
 
@@ -62 +63 @@
                 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);
@@ -72 +73 @@
                 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);
 
@@ -108 +109 @@
 
         /*clean-up*/
-        if(conserve_loads) delete loads;
+        if(conserve_loads){
+                delete femmodel->loads;
+                femmodel->loads=savedloads;
+        }
         xdelete(&uf);
         xdelete(&ug);

issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp

@@ -63 +63 @@
 
                 /*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);
@@ -77 +77 @@
 
                 /*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);

issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp

@@ -53 +53 @@
         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);